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Abstract 

A mathematically challenging model of dynamic wetting as a process of interface formation has been, for 
,-0 the first time, fully incorporated into a numerical code based on the finite element method and applied, 

P^P as a test case, to the problem of capillary rise. The motivation for this work comes from the fact that, 

as discovered experimentally more than a decade ago, the key variable in dynamic wetting flows — the 
dynamic contact angle — depends not just on the velocity of the three-phase contact line but on the entire 
flow field/geometry. Hence, to describe this effect, it becomes necessary to use the mathematical model that 

i 1 has this dependence as its integral part. A new physical effect, termed the 'hydrodynamic resist to dynamic 

wetting', is discovered where the influence of the capillary's radius on the dynamic contact angle, and hence 
on the global flow, is computed. The capabilities of the numerical framework are then demonstrated by 
I comparing the results to experiments on the unsteady capillary rise, where excellent agreement is obtained. 

Practical recommendations on the spatial resolution required by the numerical scheme for a given set of 
non-dimensional similarity parameters are provided, and a comparison to asymptotic results available in 
limiting cases confirms that the code is converging to the correct solution. The appendix gives a user- 
friendly step-by-step guide specifying the entire implementation and allowing the reader to easily reproduce 
CZ) all presented results, including the benchmark calculations. 

. . 
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1. Introduction 

Reliable simulation of flows in which a liquid advances over a solid, known as dynamic wetting flows, 
is the key to the understanding of a whole host of natural phenomena and technological processes. In the 
technological context, the study of these flows has often been motivated by the need to optimize continuous 
coating processes that are routinely used to create thin liquid films on a product [I], for example, in 
the coating of optical fibres 0[3]. However, more recently, discrete coating, in particular inkjct printing 
of microdrops [1], has matured into a viable, and often preferable, alternative to traditional fabrication 
. — processes, e.g. in the additive manufacturing of 3D structures or the creation of P-OLED displays [HUH], and 

it is becoming a new driving force behind research into dynamic wetting phenomena. In most cases, such 
flows can be regarded as microfluidic phenomena, where a large surface-to-volume ratio brings in interfacial 
effects on the flow that are not observed at larger scales. 

Obtaining accurate information about micro and nanofluidic flows experimentally is often difficult and 
usually costly so that, consequently, a desired alternative is to have a reliable theory describing the physics 
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that is dominant for this class of flows and incorporate it into a flexible and robust computational tool 
which can quickly map the parameter space of interest to allow a specific process to be optimized. Such 
computational software could be validated against experiments at scales and geometries easily accessible to 
accurate measurement and then used to make predictions in processes inaccessible to experimental analysis. 

The discovery that no solution exists for the moving contact-line problem in the framework of standard 
fluid mechanics [7J [5] prompted a number of remedies to be proposed, which are summarized, for example, 
in [ch. 3 of Ej- Of these, most are what we shall refer to as 'conventional' or 'slip' models, in which the 
no-slip condition on the solid surface is relaxed to allow a solution to exist, with the Navier-slip condition 
|10j being the most popular choice. As a boundary condition on the free-surface shape at the contact line, 
one has to specify the contact angle formed between the free surface and solid. In conventional models, this 
angle is prescribed as a heuristic or empirical function of the contact-line speed and material parameters of 
the system, e.g. (llU12j . Such models provide predictions that adequately describe experiments at relatively 
large scales, often around the length scale of millimetres. It is well established that on such scales many 
of the proposed models work equally well and that finding significant deviations between their predictions, 
and hence ascertaining which best captures the key physical mechanisms of dynamic wetting, is practically 
impossible [HCS]. 

A physical phenomenon that gives an opportunity to distinguish between different models came to be 
known under a 'technological' name of the 'hydrodynamic assist of dynamic wetting' (henceforth 'hydrody- 
namic assist' or simply 'assist'). The essence of this effect, first observed in high accuracy experiments on the 
curtain coating process [151 116) . is that for a given liquid spreading over a given solid at a fixed contact-line 
speed, the dynamic contact angle can still be manipulated by altering the flow field/geometry, for example, 
in the case of curtain coating, by changing the flow rate or the height from which the curtain falls. This 
effect has profound technological implications as it allows the process to be optimized by off-setting the 
increase of the contact angle with increasing contact-line speed by manipulating the flow conditions and 
hence postponing air entrainment |15j . 

The dependence of the dynamic contact angle on the flow field has also been reported in the imbibition 
of liquid into capillaries [T7J [TH] , in the spreading of impacted drops over solid substrates [THJ HO] and in the 
coating of fibres However, in many of these flows it is yet unclear whether hydrodynamic assist actually 
occurs, or whether the free surface bends significantly beneath the spatial resolution of the experiments, 
whereas for curtain coating the hope of attributing assist to the poor spatial resolution of experiments 
has been quashed by careful finite element simulations which show that the degree of free-surface bending 
under the reported resolution of the measurements is too small to account for the observed effect and that 
conventional models cannot in principle describe this important physical effect |21j . 

Currently, the only model known to be able to even qualitatively describe assist [22j [23] is the model 
of dynamic wetting as an interface formation process, first introduced in [24) and discussed in detail in [9]. 
This model is based on the simple physical idea that dynamic wetting, as the very name suggests, is the 
process in which a fresh liquid-solid interface (a newly 'wetted' solid surface) forms. Qualitatively, the origin 
of the hydrodynamic assist is that the global flow influences the dynamics of the relaxation of the forming 
liquid-solid interface towards its equilibrium state and hence the value of this interface's surface tension 
at the contact line, which, together with the surface tension on the free surface, determines the value of 
the dynamic contact angle. When there is a separation of scales between the interface formation process 
and the global flow, the 'moving contact-line problem' can be considered locally and asymptotic analysis 
provides a speed-angle relationship which is seen to describe experiments just as well as formulae proposed 
in other models [S]- However, in the situation where the scale of the global flow and that of the interface 
formation are no longer separated, the influence of the flow field on the dynamic contact angle will occur and 
hence no unique speed-angle relationship will be able to describe experiments. Then, the interface formation 
model becomes the only modelling tool, and, given that the processes of practical interest are free-surface 
flows under the influence of, at least, viscosity, capillarity and inertia, it is inevitable that, to describe such 
flows, one needs computer simulation, i.e. the development of accurate CFD codes, which, for the effect of 
hydrodynamic assist to be captured, have to incorporate the interface formation model. 

Since the interface formation model came into wider circulation, there has been considerable interest, 
e.g. [251 1191 120] . in using it in its entirety as a practical tool for describing dynamic wetting phenomena, 
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especially on the microscale. Review articles have also referred to the description of assist using this model 
as one of the main challenges in the field [26] . The major obstacle in the development of this computational 
tool is the mathematical complexity of the interface formation model, as one has to solve numerically the 
Navier-Stokes equations describing the bulk flow subject to boundary conditions which are themselves partial 
differential equations along the interfaces and in their turn have to satisfy certain boundary conditions at 
contact lines confining the interfaces. These conditions determine the dynamic contact angle and hence 
influence the free-surface shape, which exerts its influence back on the bulk flow. Thus, the bulk flow, 
the distribution of the surface parameters along the interfaces and the dynamic contact angle that these 
distributions 'negotiate' become interdependent, with the dynamic contact angle being an outcome of the 
solution as opposed to conventional models where it is an input. This interdependency is, on the one hand, 
the physical essence of the experimentally observed effect of hydrodynamic assist to be described but, on 
the other hand, it is this very interdependency that, coupled with the nonlinearity of the bulk equations and 
the flow geometry, is the reason why the model is difficult to implement robustly into a numerical code. 

Some numerical progress has been made for the computationally less complex steady Stokes flows [23], 
but what is lacking is a step-by-step guide to the implementation of this model in the general case, with 
unsteady effects in the bulk and interfaces as well as nonlinearity of the bulk equations fully implemented. 
This which would pave the way for incorporating the interface formation model into existing codes as well 
as developing new ones. Therefore, the first objective of the present paper is to address these issues and 
provide a digestible guide to the model's implementation into CFD codes. Then, after giving benchmark 
computational results to verify this implementation, we consider a problem of imbibition into a capillary, 
compare the outcome with experiments and predict essentially new physical effects. 

A starting point in the development of the aforementioned CFD code is to first develop an accurate 
computational approach for the simulation of dynamic wetting flows using the mathematically less complex 
conventional models and this was achieved in [27] . It was shown that many of the previous numerical results 
obtained for dynamic wetting processes are unreliable as they contain uncontrolled errors caused by not 
resolving all the scales in the conventional model, most notably the dynamics of slip and the curvature of 
the free surface near the contact line. Benchmark calculations in |27j for a range of mesh resolutions resolved 
previous misunderstandings about how to impose the dynamic contact angle and showed that implementing 
it using the usual finite element ideology, as opposed to 'strong' implementation of the contact angle, works 
most efficiently: it allows errors to be seen, and hence controlled, as the computed contact angle varies 
from its imposed value, instead of them being masked elsewhere in the code. Furthermore, in [28] . we have 
shown that numerical artifacts which occur at obtuse contact angles are present in both commercial software 
and in academic codes where, misleadingly, they have previously been interpreted as physical effects. By 
comparing computational results to analytic near field asymptotic ones, we have shown that the previously 
obtained spikes in the distributions of pressure observed near the contact angle are completely spurious, 
and, to rectify this issue a special method, based on removing the 'hidden' eigensolutions in the problem 
prior to computation, has been developed [35]. In J5S], we showed that our code is capable of simulating 
unsteady high deformation flows by comparing to benchmark calculations published in the literature and 
performed by various research groups. In contrast, in [30 it has been shown that when commercial software 
is used to simulate relatively simple benchmark free surface flows, the converged solution is not the correct 
one. 

In this series of articles [37] [551 US] j our approach has been to carefully develop a robust CFD algorithm 
for the simulation of dynamic wetting flows. Thus far, this has been achieved for the conventional models, 
where we validated our results by performing mesh independence tests, comparing with asymptotic solutions 
in limiting cases and, where no analytic progress was possible, to previously obtained benchmark solutions 
published in the literature. In this article, we continue this series of papers and, for the first time, incorporate 
the interface formation model into our code in full. Notably, it will be shown that a code implementing 
the mathematically complex interface formation equations can easily recover the much simpler conventional 
models by setting certain parameters to zero. This allows the same code to be used to compare the predictions 
of dynamic wetting models for a range of flows and hence to determine, by a comparison with experiment, 
where the physics of interface formation manifests itself in a significant way and where mathematically 
simpler conventional models are adequate. Here, we focus on dynamic wetting; the extension to other 
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flows of interest where interface formation or disappearance also occurs is a straightforward procedure 
computationally, and it will be dealt with in forthcoming articles. 

The layout of the present article is as follows. First, in f|2] without limiting ourself to a particular flow 
configuration, the equations describing the dynamic wetting process are briefly recapitulated. Then, in <|3j 
the finite element equations are derived for the dynamic wetting flow as an interface formation process. Local 
element matrices, and additional details about the finite element procedure are provided in the Appendix, 
which complements a 'user-friendly' step-by step guide to finite element simulation given for this class of 
flows in [27] and allows one to reproduce the benchmark simulations in fj4] These simulations are performed 
for the dynamic wetting flow through a capillary both in the case where the meniscus motion is steady and for 
the unsteady imbibition of a liquid into an initially dry capillary. Computations are checked for convergence 
both as the mesh is refined and towards asymptotic results in limiting cases. Having established the accuracy 
of our approach, in ^5] new physical effects are discovered by considering the influence of capillary geometry 
on the dynamic wetting process. Finally, the computational tool's ability to easily describe experimental 
data is shown in ^6] and a number of advantages over current approaches, particularly in the initial stages of 
a meniscus' motion into a capillary tube, are highlighted. Concluding remarks and areas for future research 
are discussed in Sj7] 



2. Modelling dynamic wetting flows as an interface formation process 



Consider the spreading of a Newtonian liquid, with constant density p and viscosity p,, over a chemically 
homogeneous smooth solid surface. The liquid is surrounded by a gas which, for simplicity, is assumed to 
be inviscid and dynamically passive, of a constant pressure p g . Let the flow be characterized by scales for 
length L, velocity U, time L/U, pressure pU/L and external body force F . In the dimensionless form, the 
continuity and momentum balance equations are then given by 



0, 



Re 



du 

~dt 



u • Vu 



V • P + St F, 



(1) 



where 



-pi + 



Vu+ (Vu) 



(2) 



is the stress tensor, t is time, u and p are the liquid's velocity and pressure, F is the external force density 
and I is the metric tensor of the coordinate system. The non-dimensional parameters are the Reynolds 
number Re — pUL/p and the Stokes number St — pF^L 2 /(pU). 

Boundary conditions to the bulk equations are required at the liquid-gas free surface x = Xi(si, Sz,t), 
whose position is to be found as part of the solution, and at the liquid-solid interface x = X2(si, S2, i), whose 
position is known, and at other bounding surfaces which will be specified by the problem of interest; here, 
(si,S2) are the coordinates that parameterize the surfaces. Boundary conditions along the free surface, the 
liquid-solid interface and the contact line at which they intersect are provided by the interface formation 
model [9], as follows. 



2.1. The interface formation model 

To represent the boundary conditions on an interface with a normal n in an invariant form, it is convenient 
to introduce a (symmetric) tensor I — nn, which is essentially a metric tensor on the surface. If an arbitrary 
vector a is decomposed into a scalar normal component a n — a • n and a vector tangential part a||, so that 
a = a|| + a n ii, we can see that, because n (I — nn) = 0, the tensor (I — nn) extracts the component of a 
vector which is tangential to the surface, i.e. a • (I — nn) = an. Hereafter, n is the unit normal to a surface 
pointing into the liquid, and subscripts 1 and 2 refer to the free surface and solid surface, respectively. 

The equations of interface formation on a liquid-gas free surface, which act as boundary conditions for 
the bulk equations ([I]), are given by 

(3) 
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Vu+(Vu) • (I - nini) + Vcri = 0, 

p + xix- Vu + (Vu) T ■ ni| = (TiV • ni, 
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(4) 
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(6) 
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whilst at liquid-solid interfaces formed on a solid which moves with velocity U, the equations of interface 
formation have the form 

(v| - U) ■ n 2 = 0, (10) 



Ca n 2 • P • (I- n 2 n 2 ) + ±Vct 2 = (3 (u N - U,,) 



V 2|| - 2 l U ll " 

(u - vf) • n 2 

9p| 



U, 



aV(T 2 , 



at 



V-(p 2 v|) 



= - (P2 - Pie 



C r 2 = A(l-p|). 



(11) 
(12) 
(13) 

(14) 

(15) 



The differential term <7iV • ni, where <j\ is the (dimensionless) surface tension on the free surface, 
describing the capillary pressure in the normal-stress equation ^ indicates that this equation requires its 
own boundary condition where the free surface terminates, i.e. at the contact line. There n 2 is known and ni 
must be specified by setting the dynamic contact angle 8^ at which the free surface meets the solid surface: 

ni • n 2 = - cos6> rf . (16) 

This angle is determined from a force balance at the contact line, given by Young's equation |31) : 

cr 2 + (7! cos6» d = er s , (17) 

where cr 2 and as are the surface tensions of the liquid-solid interface and solid-gas interface, respectively, 
and the latter is henceforth assumed to be negligible as ~ 0. Equations and ( 14) also require a boundary 



condition at the contact line where the two interfaces meet, and here we have the continuity of surface mass 
flux 

p\ (vJ N - U c ) ■ mi + p\ (v||| - U c ) ■ m 2 = (18) 

where U c is the (dimensionless) velocity of the advancing contact line and is the inward unit vector 
tangential to surface i and normal to the contact line (see Figure [lj . 

The interface formation model is described in detail in the monograph [S] and a series of preceding papers 
[e.g-IMlES so that here only the main ideas are briefly recapitulated. The surface variables are in the 'surface 
phase', i.e. physically in a microscopic layer of liquid adjacent to the surface which is subject to intermolecular 
forces from two bulk phases. In the continuum limit, this microscopic layer becomes a mathematical surface 
of zero thickness with p s denoting its surface density (mass per unit area) and v s the velocity with which it is 
transported. The following non-dimensional parameters have been introduced a — aa/(UL), j3 — (f3UL)/a, 
e = Ut/L, A = 7P(q)/c and Q = pf Q %/(pZ7r) which are based on phenomenological material constants a, 
(3, 7, r and P( y, in the simplest variant of the theory, the latter are assumed to be constant and take the 
same value on all interfaces; a is the characteristic surface tension for which it is convenient to take the 
equilibrium surface tension on the liquid-gas interface. 
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Figure 1: Illustration showing the vectors associated with the liquid-gas free surface A\ and the liquid-solid surface A 2 in the 
vicinity of the contact line C c ;. 



Estimates for the material constants have been obtained by comparing the theory to experiments in 
dynamic wetting, e.g. in [33] . but could just as easily have been taken from an entirely different process in 
which interface formation is key to describing the dynamics of the flow 34, 35, 36, 37 . In other words, once 
obtained for a particular liquid in one set of experiments, the material constants determined can be used to 
describe all phenomena involving the same fluid in which interface formation dynamics 'kicks in'. 

The surface tension Oi is considered as a dynamic quantity related to the surface density p\ via the 
equations of state in the 'surface phase' (|9| and (15 1, which are taken here in their simplest linear form. 
Gradients in surface tension influence the flow, firstly, via the stress boundary conditions Q, (J5j) and (11), 
i.e. via the Marangoni effect, and, secondly, in the Darcy-type equation^ Q and ( 12 ) by forcing the surface 
velocity to deviate from that generated in the surface phase by the outer flow. Mass exchange between the 
bulk and surface phases, caused by the possible deviation of the surface density from its equilibrium value 
pi, is accounted for in the boundary conditions for the normal component of bulk velocity ^ and ( 13 ), and 
in the corresponding surface mass balance equations Q and (14 1. 

One would expect a generalized set of boundary conditions to have the classical conditions as their limiting 
case. For the interface formation model this limiting case follows from the double limit e — > 0, /3/Ca —> 00. 
When the limit e — > is applied to ([3])-([9]), the liquid-gas interface equations are reduced to their classical 
form. Notably, if we apply e — > to the liquid-solid equations ( 10 )— ( 15 ) , the conventional 'slip' model is 
obtained, that is the Navier-slip condition combined with impermeability. In this limit, the surfaces are 
in equilibrium so Young's equation (17) gives that the dynamic contact angle is fixed as a constant at its 
equilibrium value 8^ = 6 e . If we wish to go further with the conventional model approach and describe the 
dynamic contact angle as some function of the contact-line speed, then Young's equation must be discarded 
in favour of this function. Therefore, implementing the interface formation model into a numerical code 
allows one to test all conventional models of wetting proposed in the framework of continuum mechanics. 
By applying the limit /3/Ca (= j3L/ 'a) ->oowe recover the classical equations of no-slip and impermeability 
on the solid surface for which the moving contact-line problem is known to have no solution [7J [5] . 

Despite the model's complexity, in limiting cases, analytic progress on it can be made to obtain explicit 
relationships for the surface variables and, with some further assumptions, even a formula relating the 
contact-line speed to the dynamic contact angle. Such formulae are a useful test of our numerical solutions, 
and we will briefly recapitulate their outcome. 



3 The analogy with the Darcy equation is that the tangential surface velocity vf, is the average velocity of the interfacial 
layer and its deviation from that generated by the outer flow is proportional to the gradient of surface tension, which is the 
negative gradient of surface pressure as p s = —a. 
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2.2. Asymptotic formulae in a limiting case 

When the contact-line motion can be analyzed as a local problem, as opposed to cases where the interface 
formation and the bulk flow scales are not separated so that manipulating the global flow influences the 
relaxation process along the interfaces, asymptotic progress is possible. A full derivation of the results we 
use may be found in [9] and references therein; here we shall just outline the main assumptions and results. 

If in the steady propagation of a liquid-gas free surface over a solid substrate in the Stokes regime 
(Re <C 1), the characteristic length scale of the interface formation process I = Ut is small compared to the 
bulk length scale L, we have that our non-dimensional parameter e <C 1. If in the limit e — > we also assume 
that the capillary number Ca — » 0, then, to leading order in Ca, the normal-stress boundary condition ([5| 
gives that the free surface near the contact line is planar, so that the problem may be considered locally in 
a wedge-shaped domain. Then, we can identify the following three asymptotic regions: 

(a) The outer region, where, in a reference frame moving with the contact line, one has a flow in a wedge 
in the classical formulation, with a zero tangential-stress and a no-slip boundary, described in |38j ; 

(b) The intermediate region with the characteristic (dimensionless) length scale e, where the surface- 
tension-relaxation process takes place and where, due to smallness of Ca, to leading order the influence 
of the bulk flow on this process can be neglected; 

(c) The inner region, with the characteristic length scale eCa, through which the surface densities and 
velocities, to leading order, stay constant. 

On the free surface, at leading order in the intermediate and inner regions, one has 



Pi = Pie 



'111 



u f (e d ), 



(19) 



where Uf(9 d ) is the (dimensionless) radial velocity of the bulk flow in the far field on the liquid-gas interface 
given by [35]: 

sin6» d - 9 d cos9 d 

sm e d cos e d - e d 

so that the surface mass flux into the contact line is —p s le Uf(9 d ). Then, since the surface variables are 



constant through the inner region, the boundary conditions at the contact line (17), (18) can be applied 



to the distributions of the surface variables in the intermediate region. As a result, we have two first-order 
ODEs to solve for the distributions of p 2 and v^n along the liquid-solid interface 



d(P2^||) 



ds 



= -(pI-pL 



(S>0) 



(21) 



where s = s/e is the intermediate region's variable and V 2 = /3e/((l+4a/3)A) is the non-dimensional contact- 
line speed, subject to (a) boundary condition (18), now taking the form p 2 v 2 = —p s le Uf(9 d ) at s = 0, (b) the 



matching condition p 2 — ¥ p 2e as s — > oo, and, as boundary condition (a) implicitly depends on the parameter 
8 d , (c) Young's equation (17), now taking the form cos9 d — X(p 2 — 1). The equilibrium contact angle 6 e is 
obviously related to p 2e by cos# e = X(p\ e — 1), which can be used to replace p s 2e with 6 e . 

The above problem is easily solved numerically; however, by taking an additional assumption A > 1, 
one can obtain an analytic relation between the contact angle and the non-dimensional speed of the contact 
line V: 



2V 



cos 9 a 



cos 9, 



+ (i - pLr 1 (i 



PleUf{9 d )) 



where 



2V(p s 2 



(V 2 + P s 2e ) 



V +[V 2 + l + cos6» e (l- pf e )] 1/2 

2V(p s 2e + pt e u f (9 d )) 



(22) 



1/2 



V 



c 



(v 2 + P s 2e ) 1/2 + v 



The presence of Uj 



in ( 22 ) highlights a connection between the flow in the outer asymptotic region 



and the value of the dynamic contact angle. When the contact line is 'insulated' from the global flow by 



7 



the low-Reynolds-number region, the flow near the contact line is completely determined by the contact-line 
speed, and hence the mass flux into the contact line that 'feeds' the liquid-solid interface can be found from 



the local solution. This is the case considered in 24, 39] , where the theory shows excellent agreement of (22 1 
with experiments. If the outer flow is to influence the mass flux into the contact line [53], this will affect the 
value of the contact angle. 

The most important length scale L,/ associated with the interface formation process is the characteristic 
distance over which the solid surface returns to equilibrium and the asymptotic result indicates that this is 

~2V(y/W+pf e -V) . 



given, in non-dimensional units, by kLif/e = 1 in (191 so that Lif 
The expressions given above will be used in Sj4] 



ow to compare our computations to in the situations 



where the underlying assumptions of the asymptotics are satisfied. Now, we shall consider the development 
of this computational algorithm for the general case without making any simplifying assumptions. 



3. Finite element procedure 

A finite element framework for the simulation of dynamic wetting flows using the conventional models 
of dynamic wetting was described in |27j . To handle the evolution of the free surface this framework uses 
an arbitrary Lagrangian Eulerian (ALE) scheme based on the method of spines, a computational approach 
which has been successfully applied to a range of coating flows over the past thirty years, e.g. in [4"01 HT1 [2Tj . 
In |2S], the framework was extended for the simulation of time-dependent free surface flows, with the code 
providing accurate solutions for the benchmark test case [42], [43] of a freely oscillating liquid drop. This 
confirmed that the implementation of the viscous, inertial and capillarity effects is accurate, even when the 
mesh undergoes O(l) deformations. 

What follows is the implementation of the interface formation equations into our framework. For a more 
detailed description of the basic components of the framework and a user-friendly step-by-step guide to 
implementing dynamic wetting flows into the finite element method, the reader is referred to (27j and in 
particular to the Appendix which makes it possible for the interested user to reproduce the results presented 
there. The Appendix in the present article provides additional details of the implementation and, in this 
sense, complements the Appendix in [57], allowing one to reproduce the results of fj4]and <|5j 

3.0.1. Problem formulation in the ALE scheme 

Consider how the equations of f|2] written in Eulerian coordinates x, are formulated for an ALE system 
where the flow domain x = X ( x > *) deforms in time. This deformation must be accounted for in the 
temporal derivatives of variables whose position in the Eulerian system is evolutionary, in particular, in the 
Navier-Stokes equations where the material derivative D/Dt transforms as 



Du du 

~Dt = ~dt 



_ du 



+ (u-c)-Vu. (23) 

X 



. . 9x 
Here, c (x, t) = ^ 



is the velocity of the ALE coordinates with respect to the fixed reference frame. It 

X 

can be seen that, as should be expected, for c = u we have a Lagrangian scheme whereas for c = the 
Eulerian system is recovered. 

Before considering temporal derivatives occurring in the interface formation equations, it is convenient to 
introduce the surface gradient V s , which is the projection of the usual gradient operator V onto the surface 
V' 5 = (I — nn) -V. An arbitrary surface vector a s is written in terms of components normal and tangential to 
the surface as a s = + a^n, where = a s • n, so that for its divergence one has V • & s = V s • + a s n V s • n. 
In particular, as described in [44] , points on the surfaces move with the normal surface velocity u*, i.e. 
according to the kinematic equation |3j), and an arbitrary tangential component which depends on the 
choice of mesh design. Then, on a given surface 



Therefore, in the ALE framework the left-hand side of the surface mass balance equations ([8]) and (14), 
become 

dp 3 -, 



dt 



v-(/>X) 



dt 



c 7 .v P7 + v-(p 7 v 7 ) 



where 7 = 1,2 refers to the liquid-gas and liquid-solid interface, respectively, and X s , = X s , ( x , t) are the 
corresponding coordinates. Then, ^ and ( 14 ) take the form 



dt 



• VX + VS ' ( P7 V 7ll ) + P'>7nV S ■ n 



7 rje 



0, (7 = 1,2), 



(24) 



where we have used that n 7 • V s p 7 — 0. Equations (24) can be rearranged to obtain 

V s 



dp: 



P 7 (v 7 ||-c 



7ll 



+ «„V S • n 7 } + (t - P° = 0, (7 = 1,2). 



In the limiting case, where the surface moves only normal to itself, so that c*m = 0, the usual Eulerian 
equations are recovered: 



dt 



V s - p T v 



4 V 7ll 



^7 r^7e 



9* 



^7 /^7e 



0, (7=1,2), 



whilst if the surface moves in a Lagrangian way, c 7 = vl, then the term under the divergence becomes 
identically zero, and we have 



dp 3 , 

— -r^v ■ v 7 || -ry y v yn * 
Op 



p 7 V 5 -v 7|| +^ 7 „V s -n. 



"7 ^76 



^+p 7 V-v 7 ) +^-# e =0, (7 = 1,2). 



(25) 



Here, as should be expected, there is no term representing the convection of surface density by the surface 
velocity, i.e. a term of the form v s • Vp s does not appear. 

Having reformulated the equations for the ALE system and introduced surface operators, we can now 
derive the appropriate FEM residuals. 



3.0.2. Forming the finite element residuals 

The defining feature of the FEM is that the computational domain V is tessellated into a finite number 
of non-overlapping elements, each containing a fixed number of nodes at which the functions' values are 
determined. Between these nodes the functions are approximated using interpolating functions whose func- 
tional dependence on position is chosen. In what follows, N p is the total number of nodes in V at which the 
pressure is determined, N M is the number of nodes at which the velocity components are to be found, N% is 
the number of nodes on the free surface A\, N2 the number of nodes on the solid surface A2, and N c the 
number of nodes along the contact line. 

The procedure of generating the finite element equations is well known and a detailed explanation of 
how this is achieved for dynamic wetting flows described by the conventional model is given in [37j, so that 
here we just give the main details. Functions are approximated as a linear combination of interpolating 
functions each weighted by the corresponding nodal value. In particular, we use mixed interpolation so that 
the Ladyzhenskaya-Babuska-Brezzi [45] condition is satisfiecT] with linear basis functions ip- s to represent 



4 Equal-order methods may also be used with stabilization, e.g. 46 47j, but these issues lie beyond the scope of this paper. 
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pressure and quadratic ones <j)\ for velocity: 



N p 

E^j( x )' 
j=l 



E 1 

j=i 



(x G V). 



In the Galcrkin finite clement method, the basis functions which are used to discretize the functions of 
the problem are also used as weighting functions to create the weak form of the problem, see [§3 of H7J for 
specific details. Note that here we use Roman letters for the indices (i, j, etc) to refer to the nodal values and 
approximating functions spanning the whole domain (globally); in the Appendix, where all the numerical 
details are given, these indices will be used alongside italicized ones etc), which will refer to local, 

element-based, values and functions. 

Surface variables are also approximated quadratically with basis functions on surface A 1 (7 = 1,2) 
denoted by 7) j(si, 82)- So, all of the interface formation variables, represented by an arbitrary surface 
variable af,, as well as the shape of the free surface Xi, are approximated as 



JV 7 

E 

.1=1 



a 7,i ( ? i '7j( S l' S 2), 



Xl 



E 

j=l 



x l,j</>l,j(si, s 2 ). 



To determine the free surface shape, i.e. the nodal values Xi.j, a function h — h(si, s 2 ) is found as part of 
the solution at free-surface every node, so that hj = hj(si,S2) for j = 1, ...,Ni, which points in a direction 
linearly independent from both si and s 2 at each node, i.e. 'out' of the free surface. For example, in the 
simplest case of a Cartesian coordinate system one could have the free surface at (x, y, z) ~ (s 1; s 2 , M s i> s 2))> 
so that h is the height above the (x, y)-plane, or, in a two-dimensional example, one may have a polar 
coordinate system with the free surface given by (0,r) ~ (si, h(si)), in which case h is the distance of the 
free surface from the origin for every angle 9. 

The basis functions used to approximated the variables are now used to derive the weak form of the 
problem, i.e. the finite element equations. From ([!]), the continuity of mass (incompressibility of the fluid) 
residuals Rp are 



R 



^V-udV (i = l,...,N p ). 



(26) 



After projecting the mom entum equations ([I]) onto the unit basis vectors of the coordinate system 
e a (a = 1, 2,3) and using (23), our momentum residuals R i ' a take the form 



R M, a 



R( |^ + (u-c)-Vu) V P ,S7 F 



dV (i = l,..,iV u ). 



(27) 



Integrating by parts and using the divergence theorem, as shown in |27j . one can rewrite (27) in terms of 
volume and surface contributions: 



+(ijf' a ) (i = l,...,N u ), 



Re 



du 



(u - c) • Vu - St F 



hi e a ■ P • n dA. 



V(0ie Q ) :P^ dV, 



(28) 



In (28), only when node i is on the surface A will be non-zero, i.e. it is nodes on the surface of V 



which contribute to the momentum residuals via ( 28 ) . This term allows us to incorporate stress boundary 



conditions naturally, i.e. by adding them as contributions to the momentum equations at the surfaces, 
which is a well known advantage of the finite element method over other numerical approaches where special 
boundary approximations have to be constructed. 
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To incorporate our free-surface stress boundary conditions into (281, equations Q and ^ are rewritten 
into the computationally favourable form 

Ca {mp g + m • P) + V s • [ax (I - rum)] = 0, 

where a\ (I — is the surface stress, playing the same role on the surface as P does in the bulk. Then 

( 28 1 can be rewritten on the free surface as 

/ 0i,ie a • P • n x dA x = / </>i.ie Q • {V s • [en (I - nim.)] + Cap g n x } dA x . 

J A! ° a J A x 

It was initially suggested by Ruschak [48], and generalized for three-dimensional problems in [49j[50], that, 
by using the surface divergence theorem, one could lower the highest derivatives and thus both reduce the 
constraints on the differentiability of the interpolating functions which are used to approximate the free 
surface shape and give a natural way to impose boundary conditions on the shape of the surface where it 
meets other boundaries. This is achieved by, first, using the chain rule: 

/ <f>i ti e a • P • ni dA x = -yr [ {V s • [<£i,iO-ie a • (I - n x ni)] - ctiV s • (<£i,ie a ) + Cap g <p 1A (e a ■ ni)} dA 1} 
JA X JA X 

and then using the surface divergence theorem |51l p. 224], which for a surface vector a s with no normal 
component, a s = af,, is given by 



V s • af, dA = - J af, • m dC, (29) 



where the unit vector m is the inwardly facing normal to the contour C that confines A (Figure [T]), with 
a fl = (fri^iea ■ (I — nini), so that 

/ 0i,ie Q P • ni dAi = — / [criV 5 • (0i,ie Q ) - Cap g (f> lti (e a ■ m)] dA x + — / 0i.iCTie Q • mi dC x . 

J Ai <^ a JAi Ca JCi 



Thus, on the free surface, the term (28) in the momentum residual is now replaced by a different surface 
contribution and a line contribution 

(af 1 ' ) = n I [fTiV s • (0i,ie Q ) - Capg<j) i;i (e a ■ m)] dA 1: (llf La ) =11 0i,iO"ie a -ini dC x . 
v i Ai ua j a x \ / Ci Ca J Cl 

(30) 

The same procedure of integrating by parts and using the divergence theorem has been used on both 
the surface stress term V s • [a\ (I — nini)] and the bulk stress term V • P. In both cases, this has created 
contributions from the boundary of that term's domain, i.e. the confining contour and surface, respectively. 
Consequently, the momentum residual now contains a cascade of scales 

R^ a = (i?, M -) v + (flf-) a + (af-) c , (3i) 

which represent, respectively, the volume, surface and line contributions. In particular, part of the contour 
C\ which bounds the free surface is the contact line C c i where the free surface meets the solid. Other 
boundaries to the free surface further away from the contact line, for example an axis or plane of symmetry, 
are treated in the same way but are not to be formulated until specific problems are considered in Sj4j 

At the contact line, it is useful to rearrange the term in the integrand of the contour integral in (30) 
by representing the vector mi in terms of a linear combination of its components parallel to ri2 and rri2 
(Figure [T]): 

m, = (mj • m 2 ) m 2 + (mj • n 2 ) n 2 . 



11 



This identity can be used to make the contribution to (31 1 coming from the contact line dependent on the 



known vector ri2 , the vector m.2 , varying along the contact line but independent of the free-surface shape, 



and 6 c i, defined by (16) and determined by Young's equation (17 1, so that 



fi? M ' Q ) = -pr I ;e Q ■ (m 2 cos fld + n 2 sin6» d ) dC t 

V 1 / c cl Ca J c 



l (i = l,...,iV c ). 



If the contact line's tangent vector is t c , then m.2 = ±t c x n2, with the sign chosen to ensure the inward 



facing vector m 2 is selected. Thus, equation ( 16 ), which defines the contact angle, can be applied in a natural 



way, that is without needing to drop another equation in order to make room for an equation that would 
fix the shape of the free surface at the contact line. The contact angle itself is determined from Young's 



equation (17), which, when put in residual form as an integral over the contact line contour, gives 



ii i (oi cos6> d + cr 2 ) dC c i (i= 1,...,JV C ). 



At the liquid-solid interface the approach developed in |27j is used. Instead of dropping the momen- 
tum equation normal to the solid to impose a Dirichlet condition on the normal velocity (13), we use it to 



determine the normal stress acting on the liquid-solid interface which allows the contact line, where bound- 
ary conditions of different type meet, to be treated in a manner consistent with standard FEM ideology. 
Specifically, we introduce a new unknown A r\ on the liquid-solid interface which is defined by the equation 



A = n 2 • P • n 2 . 



(32) 



It should be pointed out that this particular implementation simplifies the finite element procedure indepen- 
dently of the dynamic wetting model chosen and is particularly useful when considering a surface non-aligned 
with a coordinate axis. In this case, the procedure of rotating the momentum equations to align with the 
coordinate axes [52] is cumbersome whereas our approach is independent of both the free surface and the 
solid's shape, that is the curved nature of a surface is as easy to handle as, say, a planar surface aligned 
with coordinate axes. 

On the liquid-solid interface the contribution to the momentum equations from the stress on the surface, 



which contains contributions from both the generalized Navier condition (11) for the tangential stress and 



( 32 ) for the normal stress gives 



>2a 



A (e Q • n 2 ) + — e Q • (u M - U||) - 7^ e a ■ W 2 



2Ca 



dA 2 



(i=l,...,N 2 ). (33) 



where 02, i is an interpolating function for the liquid-solid interface corresponding to the i-th node. 

In addition to the boundary conditions involving stress on each interface, we have an additional equation 
involving the velocity normal to each interface. On the liquid-gas free surface this is the kinematic condition 
pi) whose residuals are given by 



dAi 



(i = l,...,iVi), 



(34) 



whilst on the liquid-solid interface we have a condition of impermeability of the solid (10) with residuals R[ 
given by 

R(= [ 02,i (v|- U) • n 2 dA 2 (i = l,...,N 2 ). (35) 
Ja 2 

In keeping with the framework presented in |27j , all momentum equations are applied at both the liquid- 
gas free surface and at the liquid-solid interface and hence, once the two boundaries meet at the contact line, 



5 Labelled A in [271. 
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the contact line conditions can be implemented naturally, without dropping any of the equations there. A 
crude, but useful, interpretation is to think of the momentum equations as determining the bulk velocities, 
the kinematic condition on the free surface as determining the unknowns that specify its position, and the 
impermeability condition, which is the geometric constraint of the prescribed shape, as determining the 
'extra' unknown A, i.e. the normal stress. 

Thus far, the equations are assumed to determine the bulk velocity, the shape of the free surface and 
the normal stress on the liquid-solid interface (A). In addition to these equations we have on the 

free surface and ( 12 )— ( 15 1 on the liquid-solid interface which determine the surface velocity v s , surface 
density p s and surface tension a. In particular, we can think of the Darcy-type equations ([6| and (12) as 
determining the surface velocity tangential to the surface vf,, which in a fully three-dimensional flow will 
have two components in the e Sl and e S2 direction, where e s „ is a basis vector in the direction of increasing 

s u , v = 1,2. The residuals i?; 1 " ' andi?; 2 " are 



R 11 



ill 



U|| 7jr — Wx ) dA 1 



4/3 



R? U ' V = I 4>2,ie Su ■ (v|n - | (u||+U||) -svv 2 



dA, 



(i=l,..,Aq), 



(i= 1,...,A 2 ). 



(36) 



(37) 



Equations ([7]) and ( 13 ) can be thought of as determining the normal component of the surface velocity on 
the surface A 7 (7 = 1,2) and the residuals i? ; 7 " from these equations take the same form on each interface: 



R 



7ii [(u-v;).n 7 -g(p;- P ; e )] dA, (i = l,..,iV 7 ; 7 = 1,2). (38) 

l^) and |l4| to determine the evolution of the surface 



The corresponding residuals R i 1 from equations 
density are 



R 



dt 



V 7ll Sll 



+p 7 V s ■ c 7 || + P ;(v 7 • n 7 )V s • n 7 ) + p 7 - p 7e } 



dA-, 



(i = l,...,/V 7 ; 7-1,2). (39) 



Using the standard FEM ideology, which will give us a method for applying boundary conditions on the 
surface flux p s v|j where the surface meets a boundary, we integrate the divergence term in (39) by parts to 
obtain 



& = 



-ep: 



^7,^7 ( 



7 I V 7ll ~ Sll 



tJJ 



■V s 



"7 1 1 



^7,1 



^7,1 



dp* 



p'V 3 • c* I, + p 7 (v 7 • n 7 )V s • + pt - p 



dA-, 



gt i rj - — r| I T^ 7 iv 7 • xx 7 

(i = l,...,/V 7 ; 7 = 1,2). 



ye 



gL4, 



Then, using the surface divergence theorem (29) with = 4>\ t \p s (v| — C|| 



both the surface and the contour bounding that 



l ll 

surface, so that 



(40) 

we have a contribution from 



(#1 



dp^ 

dt 



P^ s -c 



Til 



- P y in V s • n 7 + pt] (v 7 • n 7 ) V s 



dA 



e07,i-° 7 ( V 



Til Sll 



m 7 dC 1 



(1 = 1, 



•,JV 7 ; 7 = 1,2) 



(41) 

7! 

(42) 
(43) 



13 



Consequently, we are able to apply any boundary conditions on surface mass flux by replacing the contour 
contribution (431 with the appropriate value. In particular, at the contact line contour C c ; we have the 
surface mass flux continuity condition (18). This condition could potentially be applied by using it to 



replace either the contour contribution to the free surface or the liquid-solid interface equations, and to 
determine which way is correct we study the structure of the interface formation equations. Specifically, an 



asymptotic approach to the dynamic wetting problem for small Ca and e recapitulated in \ 2.2 shows that 
it is the flux of surface mass into the contact line from the free surface, which depends on the global flow 
via the velocity in the far field on the free surface Uf(0d), that determines the relaxation process on the 
solid surface. Therefore, in a numerical implementation of this problem one should allow the free surface 



equations to determine the flux into the contact line and then use the surface mass continuity condition ( 18 1 



on the liquid-solid side of the contact line to take this flux as the surface mass supply into the liquid-solid 
interface. 

Therefore, on the free surface we have 



(V) Cc; =-f ^UPl (v?,| -c?||) -mi dC c , 



whilst on the liquid-solid interface, using ( 18 ) to rewrite the flux p| ^v. 
terms of the flux into the contact line from the free surface, we have 



2|| 



-2|| 



( R f 2 ) Cci = j c e &.iP? ( V l|| - C l||) ' m l dC < 



(44) 

rri2 into this interface in 
(45) 



Then, the flux into the liquid-solid interface is given in terms of the flux that goes into the contact line from 
the free surface. 

The final residuals from the surface equations, R^" 1 , are obtained from the surface equation of state (jSJ) 
and (15), which express the surface tension as an algebraic function of the surface density: 



R 



r ,i (<r 7 - A (1 -/)$)) cL4 7 (i 



1 AT • 



7 = 1,2). 



(46) 



In this section, we have derived the finite element residuals required to solve dynamic wetting flows using 
the interface formation model. In particular we have the following coupling of unknowns 

(ui • e a ,pi,hu Ai, v 7|M • e,,, v 7nil , p 7)i , <r 7! i, d ) a = 1,2,3 7 = 1,2 v = l,2, (47) 

each of which has its corresponding residual 

(i?^^i?p,i?f,i^f,^ ll •^i^^^i^f^ J R^^fi^) " = 1,2,3 7 = 1,2 ^ = 1,2. (48) 

The subscripts i above will have different limits that are the same in both the variable and its corresponding 
residual, i.e. the same number of equations and unknowns has been assured. 



4. Validation of the code and some benchmark calculations for Problem A 

Computations are performed for cases which are axisymmetric or two-dimensional in simple geometries 
so that calculations may be easily reproduced, thus giving benchmark results for future investigators. In 
what follows, we consider a meniscus rising against gravity through a cylindrical capillary of radius R. The 
computational domain is a region in the (r, z)-plane, and the free surface is parameterized by arclength 
s. Motion is assumed to be axisymmetric about the z-axis which runs vertically through the centre of the 
capillary with the radial axis perpendicular and located at the base of the capillary (Figure [2]). First, we 
consider the steady propagation of a meniscus through a capillary (hereafter Problem A) to check convergence 
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Figure 2: Illustration showing the computational domain for flow through a capillary with the bulk domain V, liquid-gas free 
surface Ai, liquid-solid surface A2 and contact line C c ; all indicated. 

of our code as the mesh is refined and to compare our computations to the asymptotic results summarized 
in j ]2.2| Then, we will examine Problem A from the view point of analyzing the influence of the radius R 
of the capillary on the interface formation dynamics. Finally, the unsteady imbibition of a liquid into a 
capillary will be considered (hereafter Problem B) and the results compared to published experimental data 
for this type of flow. 

In the (r, z)-plane where the solution is computed, in addition to the equations formulated in Sj2j on the 
axis of symmetry for both Problem A and B we have the symmetry conditions in the form of impermeability 
and zero tangential stress 

u • n = 0, n • P • t = 0, (r = 0, < z < z a ), 

where z a is the a-priori unknown apex height. Additionally, at the apex we have the conditions of (a) 
smoothness of the free-surface shape and (b) the absence of a surface mass source or sink: 

n • e r = 0, pfv*|| • e r = 0, (r = 0, z = z a ). 

Before doing the calculations we need to consider estimates for the model's parameters for the two 
problems to be studied, leaving free the radius of the capillary as a parameter, whose influence will be 
examined in <|5j 

4-1. Typical parameter regime 

To obtain estimates for our parameters, consider the flow of a water-glycerol mixture through a capillary 
of radius R at speed U = 0.01 m s _1 . At 60% glycerol this gives fluid properties of p « 10 3 kg m~ 3 , 
p ss 10~ 2 kg m _1 s _1 and a ss 7 x 10~ 2 N m _1 O n the molecular scale, the interface is a layer of finite 
thickness I and, as discussed in |33| , the generalized Navier equation and Darcy-type equation are analogous 
to what one would have for the averaged quantities for flow in a channel of width I. Using this analogy, 
one would have /3 ~ p,/£ and a ~ £/ p, so that, taking the coefficients of proportionality as unity as a first 
approximation, and estimating p\ e ss 0.6 (so that A — 2.5), one has from dynamic wetting experiments [33] 
that the relaxation time scales as t = r M /i with extracted from the data as r M w 7 x 10~ 6 kg -1 m s 2 . 
The equilibrium contact angle is taken as 9 e = 30°, so that p| e = 1.35. Using these estimates, we have the 
following parameters of what we will regard as the base state. The parameters independent of the length 
scale are 

Ca = 10~ 2 , Q = 0.1, p s lE = 0.6, pl e = lM, A = 2.5, d e = 30°, 
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and those dependent on the length scale are 



Re=10- 2 (iim)- 1 R, St = lO' 5 (fim)- 2 R 2 , a = lO -1 ^)^ 1 , /3 = lO^m)" 1 ^, e = lO^/xinJir 1 . 

As capillary sizes of interest can range from the millimetre scale R = 10 3 fim, relevant, for example, for 
applications in microgravity |53| . right down to a few tens of nanometres R = 10~ 2 /jm [TB], the flow can be 
dominated by different physical factors. In particular, in large capillaries bulk inertia becomes important 
as Re > 1 and the surface tension relaxation effect is localized (e <C 1), being important only close to 
the contact line, whilst for nanoscale flows inertia becomes negligible (Re <C 1) and the surface tension 
relaxation length becomes comparable to the capillary's width e w 1. 

Now, we will present benchmark calculations for the steady propagation of a meniscus through a capillary 
(Problem A) and then consider the time-dependent rise of a meniscus into a capillary (Problem B), which 
will be compared to experiments. 

Jf.A.1. Problem A: steady propagation of a meniscus through a capillary 

For the steady propagation of a meniscus, the 'base' of the domain (z = 0) is sufficiently far (quantified 
below) from the meniscus so that the base can be considered as a 'far field'. In the far field, the flow is fully 
developed and the surface variables take their equilibrium values. 

The velocity distribution across the capillary base is a Poiseuille profile adapted to allow for slip at the 
liquid-solid interface that follows from the Navier-slip condition and an additional flux of mass out of the 
domain via the liquid-solid interface which occurs when Q is non-zero: 

u = 0, «, = -! + 1/2 + {2 \ + q ] mCa) (1 + Vifi/Ca) - r 2 ) , , = 2e QpL . 

The derivation of this condition is given in the Appendix. Notably, if q = and /3 jCa — > oo, then w = 1 — 2r 2 
which is the usual Poiseuille profile satisfying dw/dr = at r = 0, u> = — 1 at r = 1 and J r _ w r dr = 0. 

Taking the radius of the capillary as R — 100 //m, the parameters dependent on the capillary width 
become Re = 1, St = 10 _1 , a — 10~ 3 , /3 = 10 3 and e = 10 -4 . In this case, e = 10~ 4 is so small that the 
distance of the far field is determined by the need for the bulk flow to settle as opposed to for the surface 
phase to relax to its equilibrium state. It is known that putting the contact line a distance of ten radii away 
from the 'far field' is more than enough for this to be satisfied [Ml 12"?] . 

The streamlines for this flow are shown in Figure [3] It can be seen that the motion of the solid with 
respect to the meniscus drags fluid away from the contact line and that this in turn leads to a flux of fluid 
into the contact-line zone from the region near the free surface. In Figure |4j one can see that, as predicted 
by experiments, the motion of the liquid in the vicinity of the contact-line is 'rolling', with the liquid-solid 
interface formed by adsorbing fluid from the bulk. The Poiseuille profile, indicated by parallel streamlines 
in Figure [3j is established relatively quickly, so that the truncated far field is certainly far enough away from 
the contact-line region to not influence the results. The mesh is graded, with small elements near the contact 
line and larger ones further away, so that the essential physics of wetting, occurring on a smaller length scale 
than the bulk flow, is fully resolved whilst the problem remains computationally tractable. The details of 
the construction of this mesh are given in [27] , where it is shown how the bipolar coordinate system can be 
utilized to provide circular spines near the contact line and straight ones further away to match with the 
domain's shape. 

4-1-2. Mesh independence tests 

The mesh is constructed such that the number of elements in the body of the capillary and those near the 
meniscus can be separated. This allows us to study the convergence of the surface equations to a solution 
without being concerned by the bulk dynamics. The number of nodes along the surfaces can be chosen as 
well as how fast the element size is increased as one moves away from the contact line. Here, we fix the 
increase in element size as 10%, and check that the solution converges as the mesh is refined, i.e. as we 
increase the resolution near the contact line. 
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Figure 3: Left: streamlines for the base state simulation in increments of 0.02. Right: corresponding finite element mesh in 
this region. 




Figure 4: Left: a magnified view of the flow near the contact line for the base state simulation showing the adsorption of the 
fluid into the liquid-solid interface. The length scale is eCa = 10~ 6 and streamlines are in increments of 2.5 X 10 — s eCa. Right: 
corresponding finite clement mesh in this region. 
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0.2 0.4 0.6 0.8 1 

Figure 5: Convergence of the free-surface shape with curves corresponding to smallest element sizes of 1: 6 X 10 — 3 , 2: 3 X 10 — 4 , 
3: 2 X 10 — 5 , 4: 9 X 10 — 7 , 5: 5 X 10 — 8 . The dashed line shows the spherical cap approximation which meets the solid at contact 
angle 8 app = 72.4°. 

In |27j . it was shown that the dynamic contact angle 9d imposed into the weak formulation can differ 
from the computed one C , i.e. the one obtained from the computed free-surface shape, if the mesh is 
under resolved near the contact line. Rather than imposing the angle in the strong form, and moving the 
errors to other less observable parts of the numerical scheme, as has often been the case in previous works, 
it was demonstrated that the difference \9d — 9 C \ provides a simple error-indicator for the scheme. When 
considering the conventional model, it was showrj^] that to resolve both the dynamics of slip and the free- 
surface shape near the contact line, so that for the computed angle one has, say, \9d — c \ < 0.1°, the smallest 
clement size l m in must satisfy l m in < /3 -1 min(5 x 10~ 2 , Ca). Then, for the flow considered, this imposes 
Imin < 10 -3 min(5 x 10~ 2 , 10~ 2 ) = 10~ 5 . In addition, we now have length scales associated with the physics 
of interface formation; in particular, the asymptotics for Ca, e <C 1 shows that, in the cases considered, the 
inner most region is on the scale 0(eCa) — 10~ 6 . 

To study convergence, it is convenient to consider the following two angles: 

9d'. The dynamic contact angle featuring in the problem statement; this angle is imposed into the code 
through the weak formulation, but in the mathematical problem it is a variable whose converged value 
is to be found. 

9 C : The angle which the free surface computed for a given spatial resolution of the mesh makes with the 
solid. 

Consider 13 meshes with smallest elements ranging over five orders of magnitude from l m i n = 5x 10 
down to l m in = 5 x 10~ 8 . Coarser meshes were unable to provide a solution as the computed angle 9 C 
increased past 180°, despite the imposed angle 9d being less than 60°, in the same way as observed in [27]. 
The change in the free-surface shape during the mesh refinement procedure is shown in Figure [5j where one 
can see just how bad the approximation is on our most coarse mesh, represented by curve 1. It can be seen 
that convergence appears to take place in two stages, so that decreasing the smallest element from 3 x 1CP 4 
(curve 2) to 2 x 10~ 5 (curve 3) has little influence on the free-surface shape, but, for the reason explained 
below, these shapes are yet far from the converged free-surface shape (curve 5) which is obtained only when 
the mesh elements decrease much further. 

The two-stage convergence can be explained by examining the aforementioned angles as the mesh is 
refined, as shown in Figure § One can see that at around l m in = 10 5 the computed angle 9 C becomes 



6 Note that in this previous publication the non-dimensionalization is slightly different with /3 = flL/fi whereas here it is 
0/Ca = /3L/fi ' 
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Figure 6: Convergence of the dynamic, computed and apparent angles as the mesh, with smallest element size l m int 
for the base state simulation. The converged values of 9^ and 9 a are = 60.5° and Q a = 72.3°, respectively. 



is refined 



indistinguishable from the dynamic angle 9d that goes into the weak formulation, but at this stage 64, which 
is a part of the solution, is yet to converge to its final value. This corresponds to curves 2 and 3 in Figure [5j 
More specifically, we see that \9 d — 9 C \ < 0.1° at l min = 7 x 10~ 6 , close to 10~ 6 predicted from [27], i.e. the 
mesh is already sufficiently resolved for the computed angle 9 C to be close to the angle 0d that features in 
the mathematics of the problem and goes into the weak formulation. However, in our problem, the angle 9d 
is itself a variable and the deviation of 6d from its final value 6^ = 60.5° is still large \9d — &d\ = 3.8°. The 
convergence towards a final solution is then dominated by the requirement to resolve the interface formation 
dynamics. We have \9 d — @d\ < 0.1° at l m i n = 4 x 1CP 7 which again falls close to our estimate and shows 
the need to resolve the length scale eCa — 10~ 6 . This suggests a new estimate which builds upon that given 
in [37], that for errors in the angle of the order of 0.1°, we require 

lmin < min(5 x lO -2 /? -1 , Ca/T 1 , 10 -1 eCa), (49) 

which, respectively, ensure that the code resolves the free surface near the contact line (particularly at high 
Ca) , the length scale on which slip occurs and the dynamics of interface formation. 

A frequently used way of simplifying the problem of modelling the flow through a capillary and of easily 
interpreting theoretical and experimental results for this flow is to approximate the free-surface shape as a 
spherical cap [55] [56] . Then, the free-surface shape is fully specified by the difference between the contact- 
line height and apex height h = z c — z a . Usually, the spherical-cap approximation is characterized by the 
so-called 'apparent' contact angle 6 a , which is the angle that a spherical cap fitted through the apex and 
the contact line makes with the solid: 

a = n - arccos (- ^ ^ 2 j . (50) 

A spherical cap is the free-surface shape obtained in the limit of Ca — > and by computing 9 a we will also 
have a measure of the deviation from this shape due to viscous bending of the free surface. 

In Figure[5]the spherical cap approximation for the converged solution is plotted (dashed line). Deviations 
from the actual free-surface shape (curve 5) can be observed and the apparent angle is considerably larger, 
by 11.8°, than the dynamic contact angle, highlighting the finite capillary number regime considered, which 
will increase at higher capillary numbers. Figure [6] shows that the convergence of the value of 9 a and its 
deviation from the final value, i.e. the one corresponding to the converged free-surface shape, closely follow 
the trend of 9 C . This demonstrates that the shape of the free surface near the contact line is what governs 
the global error to the free-surface shape, as expected for small capillary numbers. In other words, the 
dynamics which governs 9d must be resolved not only to determine the flow local to the contact line, but 
also to accurately predict the entire free-surface shape and hence the global flow. 
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Figure 7: Distribution of the surface variables along the free surface against the scaled distance from the contact line s/e. 
Curves are for (Ca,e) = 1 : (5 X 10" 2 ,5 X 10" 4 ), 2 : (10" 3 ,2.5 X 10" 2 ), 3 : (10" 3 ,5 X 10" 4 ). Left: pf compared to its 
asymptotic (equilibrium) value (dashed line). Right: surface velocity tangential to the free surface and bulk velocity tangential 
to the surface compared to far-field free-surface velocity Uf = —0.68 (dashed line). 



4-1-3. Convergence to the asymptotic solution as e — > 0, Ca — > 

Having established that our computed solution converges under mesh refinement, we now check that this 
solution converges to the asymptotic results outlined in §2.2| as e,Ca —> 0. As in the the asymptotics, the 
parameter V 2 = e/3/(A(l + 4a/3)) is regarded finite in the asymptotic limit considered. To check convergence 



to the asymptotics, the parameters used in { 4.1.2 will remain unchanged except that we now vary e and Ca 
and, for a start, fix V 2 — 0.1 so that /3 = u~ x — 1.25e _1 . 

Three data sets are considered, with (Ca,e) = 1 : (5 x 10~ 2 ,5 x 10" 4 ), 2 : (10~ 3 ,2.5 x 10" 2 ), 3 : 
(10~ 3 ,5 x 10~ 4 ), so that we can ascertain the influence of decreasing either e or Ca by a factor of 50. 
The simplest 'integral' measure of convergence is the dynamic contact angle whose asymptotic value for 
this set of parameters is 102.1°. The obtained solutions for 8d for 1 — 3 are 104.1°, 107.7° and 102° so 
that the deviations from the asymptotic result are 2°, 5.6° and 0.1°, respectively. Therefore, we have the 
convergence of the numerically obtained dynamic contact angle to the asymptotic one, with an indication 
how this convergence is controlled by e and Ca. 

To examine the issue of convergence in more detail, consider the distributions of the surface parameters 
and the bulk velocity along the interfaces. Deviations of the surface density on the free surface pf from its 
equilibrium value pf e and the distributions of the tangential components of the surface velocity v{ t and the 
bulk velocity tangential to the free surface uu are shown in Figure [7J As one can see, the magnitude of the 
variation in pf along the free surface is governed by Ca, but it is e that determines how far the 'level' of pf 
is from its equilibrium value. As e — > 0, Ca — > 0, we have that pf — > p\ e . 

In agreement with the asymptotic results, the surface velocity vf t coincides with u\ t in the intermediate 
region and the two velocities begin to deviate from each other at around s/e ~ Ca, i.e. in the inner 
region, where, as expected, vf t remains constant whereas u\ t decreases as the contact line is approached. 
This pattern is clearly seen for data set 3, where both the surface and bulk velocity attain the value of 
Uf = —0.68 in the outer region, s/e 3> 1, then, as the intermediate asymptotic region is entered, vf t and 
u\t remain at this level through this region Ca <C s/e <^ 1, and as the inner region is reached, the surface 
velocity v\ t remains constant up to the contact line, whereas u\t varies. Since pf is also constant through the 
inner region, it is the surface mass flux of the intermediate region that goes into the contact line and hence 
determine the mass flux that feeds the liquid-solid interface. Notably, as one could already see in Figure |4j 
for Q 7^ the contact line is not a stagnation point for the bulk flow, so that Ui t ^ at the contact line. 

The relaxation of the surface variables along the liquid-solid interface to their equilibrium far-field values 
is shown in Figure [8] As one can see, the computed curves converge to the asymptotic distributions of §2.2| 
as e, Ca — > with data set 3 being practically indistinguishable from the asymptotic curves (dashed lines). 
In all cases the surface density is below its equilibrium value at the contact line, and it takes a finite distance 
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Figure 8: Distribution of the surface variables along the free surface against the scaled distance from the contact line s/e 
compared to the asymptotic distribution (dashed lines). Curves are for (Ca, e) = 1 : (5x 10 -2 , 5x 10~ 4 ), 2 : (10 -3 , 2.5x 10 -2 ), 3 : 
(10 — 3 ,5 X 10~ 4 ). Left: surface density Right: surface velocity tangential to the solid i>| ( . 
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Table 1: Comparison of computed solution with asymptotic value and formula ]22t. 



for the interface to fully form. At the contact line, the surface velocity tangential to the solid is lower than 
its equilibrium value and it increases to reach the solid's velocity on the same characteristic length scale. 

The convergence of our computational results to the asymptotic ones in the limit considered has now 
been confirmed. Now, we can do the reverse and use the code to ascertain the accuracy with which the 



boundary- value problem for the ODEs (21 ) gives the dynamic contact angle and the accuracy of the formula 
(22 ) obtained under the additional assumption A 1 when both are applied in a range of parameters where 
the assumptions of the asymptotics may not be satisfied. In other words, it is worth checking how the 
asymptotics works outside its formal limits of applicability. 

In Table[l] a range of values of e and Ca are considered with the parameters as above except for ft 



a 



10 3 , which is now fixed, so that V 2 varies. In general, over the considered range the asymptotics gives a 



remarkably good approximation to the dynamic contact angle. The formula ( 22 ) also gives a very reasonable 
approximation which improves considerably for smaller values of e. This good agreement is obtained despite 
the fact that V 2 is smaller than the value of Ca in certain cases, i.e. a parameter considered finite in 
the asymptotics is smaller than the one assumed to be asymptotically small, e.g. for e = 10 -4 we have 
V 2 — 8 x 10 -2 , and yet for Ca = 10 _1 the computed angle is still only 1.6° away from the asymptotic value. 



Thus, it has been shown that our computational scheme, which incorporates the interface formation 
model, is mesh-independent once all the length scales are resolved and converges to the asymptotic results 
in the limiting case. The computational results reported so far can be used as benchmarks for those who 
are interested in putting the interface formation model into a numerical code. It is useful to know that 
the asymptotic result gives a good prediction of the key feature of the wetting process, the value of the 



dynamic contact angle, with the formula ( 22 ) giving a reasonable approximation, well outside their strict 
limits of applicability. Now, we will use the developed numerical tool to make novel predictions on the 
size-dependence of the dynamic wetting flow through a capillary. 
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5. Size dependence of flow characteristics in Problem A 

Now, consider how the flow characteristics of steady capillary rise (Problem A) depend on the radius 
of the capillary. We will vary the radii over three orders of magnitude, from the R = 100 /im considered 
previously down to R = 0.1 fim. In all cases, the far field is placed at a sufficient distance from the contact 
line, ensuring that the liquid-solid interface is in equilibrium there and that moving it further away does not 
alter the dynamics of the interfacial relaxation process. 

Conventional theories, where the dynamic contact angle is an input, give that, because the material 
properties of the system remain unchanged and the meniscus is travelling at the same speed, the dynamic 
contact angle must be the same in all of these cases, i.e. independent of R. In contrast, in the interface 
formation model, where the dynamic contact angle is an output determined by the relaxation process along 
the interfaces, i.e. it is part of the solution, changes to the flow field near the contact line due to changes in 
R can alter the angle. 

In Figure [9j the bulk and surface velocities tangential to the free surface are plotted against s/e = 
10 2 (/im) _1 i? s. In the geometry of a capillary, the apex is a stagnation point in the flow; it slows down the 
flow in its vicinity so that, as one can see in the Figure, the smaller R becomes, the lower the bulk velocity is 
along the free surface. As the contact line is approached, the surface velocity follows the bulk velocity until 
s ss lOe and then the bulk velocity decreases whilst the surface velocity remains approximately constant. 
Thus, when the capillary is relatively large (R > 10/im), as for curves 1, 2, there is enough space for the bulk 
velocity and surface velocity to increase from zero at the apex to a 'far field' value, which is not far from 
u f{®d) = —0.57 given by the asymptotic solution. Then it is this value that the surface velocity keeps up to 
the contact line. However, as the radius of the capillary decreases and the apex gets closer to the contact 
line, a 'far field' no longer exists and the velocities do not have enough room to accelerate to reach the same 
'far field' value, and, as a result, the surface velocity at the contact line becomes much reduced. As can be 
seen from Figure [9j the surface density on the free surface is less affected, remaining close to its equilibrium 
value even for the smallest capillary and consequently it is the surface velocity which determines the surface 
mass flux into the contact line. For smaller capillaries this flux is significantly reduced. 

Continuity of surface mass flux across the contact line means that, as R goes down, the flux into the 
liquid-solid interface is also reduced, as can be seen in Figure [Toj where, as is the case for the free surface, 
changes in the surface velocity are more pronounced than those in surface density. The general trend for 
both of these surface variables along the liquid-solid interface is the same exponential profile, with the value 
at the contact line differing in each case. Notably, as the capillary size decreases, the influence of slip on the 
velocity at the liquid-solid interface becomes more pronounced and this accounts for curve 4 of the surface 
velocity deviating from the other curves for large s/e. 

The value of the dynamic contact angle is 'negotiated' (via Young's equation) by the surface tensions 
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Figure 10: Distribution of surface variables along the solid surface in capillaries of sizes X.R = 100 /im, 2:R = 10 fim, 
3:R = 1 /im, 4:R = 0.1 /im, against the scaled distance from the contact line s/e. Left: surface density p|. Right: velocities 
tangential to the solid surface from the bulk U2t and in the interface ti| t . 
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Figure 11: Dynamic contact angle as the capillary radius R is varied from 100 nm up to 100 fim. The lower dashed line is 
the asymptotic value using obtained from | |20[ | whilst the upper one is obtained using Uf = 0. 

at the contact line, and, although the changes in the surface density, and hence surface tension, are less 



pronounced than those in surface velocity, it can be seen in Figure 11 that the former still lead to appreciable 
changes in the dynamic contact angle as the capillary radius is reduced. There is a general trend upwards 
which is consistent with previous studies [§] showing that, if the forming liquid-solid interface is 'starved' of 
surface mass flux coming from the free surface, this results in higher dynamic contact angles. In Figure [TT] 
the lower dashed-line corresponds to the asymptotic value, whilst the upper dashed-line is the value for the 
angle which would be obtained if the velocity in the far field is taken as Uf =0 in the asymptotic formula 



(22 1. Our curve is wedged between these two potentially limiting values, and it will be the subject of further 
work to investigate whether this simple estimate is actually a bound on the variation of the dynamic contact 
angle in this geometry or whether additional factors come into play in other parameter regimes. 

6. Problem B: Unsteady imbibition into a capillary 

Next, the time-dependent imbibition of a liquid into an initially dry capillary is considered, i.e. the case 
in which a capillary comes into contact with a reservoir of liquid leading to the fluid being sucked into 
the capillary until the meniscus reaches Jurin's height, at which gravitational forces balance the capillary 
ones that drive the flow. The capillary is aligned vertical to the gravitational field so that the body force 
F = —ge z . It is assumed that the base of the capillary is maintained at a constant pressure, the atmospheric 
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one, and that the flow at the inlet is parallel to the walls of the capillary so that 



u = 0, P = P g , (0 < r < 1, z = 0). 

As an initial condition, we assume that a small amount of liquid is in the capillary, of height one-tenth of the 
radius of the capillary h — 0.1, that this liquid is at rest, and that the meniscus has not had time to adjust 
to its surroundings so that the interface is flat with z\ = h for < r\ < 1. The interface formation variables 
on the free surface are taken to be the equilibrium values, so that o\ = 1, and the liquid-solid interface is 
assumed not to be formed yet so that p s 2 = p s ^ Q y This corresponds to <72 — and hence is consistent with 



our assumption of a flat interface as the Young equation (17) gives 6 a = 90°. Thus, as initial conditions we 
have 

u = 0, Pi = Pi e , P2=P\ ), ft = 0.1. 

As a benchmark test case, parameters are obtained from the experiments in 57 1. In particular, we 
compare our simulations to experimental data in Figures 6, 8 of |57j . where the imbibition of a silicon 
oil of density p — 9.8 x 10 2 kg m -3 , viscosity p = 12.25 kg m _1 s _1 and a surface tension with air 
of a = 2.13 x 10~ 2 kg s~ 2 into capillaries with perfectly wettable walls, 8 e = 0°, was studied. The 
two capillary radii studied were R = 0.036 cm and R = 0.074 cm. Rather than fitting the interface 
formation parameters for this liquid-solid combination, we use, as a first approximation, estimates for the 
model's parameters considered previously. Using these parameters and taking as a characteristic velocity 
U = a/p = 1.74 x 10~ 3 m s" 1 , so that Ca = 1, we have for R = 0.036 cm that 

i?e = 5xl0" 5 , St = 5.8 x 10~ 2 , e = 4.1xl0" 4 , a = 5.6 x 10~ 6 , j3 = 1.8 x 10 5 , 
Q = 1.3xl0~ 2 , p s le = 0.6, p s 2e = 1.4, A = 2.5, 9 e = 0°, (51) 

whilst for R = 0.074 cm 

i?e = lxl0~ 4 , Si = 2.5 X 10" 1 , e = 2 x 10~ 4 , a = 2.7 x 10~ 6 , (3 = 3.7 x 10 5 , 

and the remaining parameters are the same. 

A relevant reference case for comparison is the so-called Lucas- Washburn approximation [5S] , often used 
to analyze experiments |58j on the unsteady imbibition into capillary tubes. In this approximation, the 
Navier-Stokes equations are volume averaged to obtain the meniscus height h as a function of time under 
the assumptions that (a) Poiseuillc flow occurs throughout the entire capillary, (b) that the driving force 
for the flow is the difference p g — p c in pressure between the pressure at the inlet p g and the pressure at the 
meniscus p Cl (c) the meniscus is a spherical cap at all times so that p g —p c = 2a cos 0d/R and, in the original 
form of the theory, (d) the contact angle is assumed to be equal to its equilibrium value the entire time, so 
that, in our case, 0d = 0° and hence p g — p c — 2a/ R. When there is no gravitational field, the meniscus 
continues to propagate indefinitely into the capillary, being slowed by the gradually increasing viscous drag 
as the distance between the meniscus and the inlet increases. If gravity is taken into account, then the 
meniscus will eventually reach its equilibrium Jurin height hoo calculated from balancing the driving force 
Pg ~ Pc with the body force on the fluid pgh aol i.e. = 2a / (pRg) . Inertial effects can also be easily 
included, resulting in a nonlinear ODE to be solved for the meniscus height which in dimensional form is 
given by 

W'+M^-lM'+^-M, (52) 
where the prime denotes differentiation with respect to time and, to compare with our calculations, the 



initial conditions at t = are h = O.li?, h! — 0. The terms on the left hand side of (52 1 are inertial forces, 
the first term on the right hand side represents viscous forces, the second one is the capillary driving force 
and the last term is the body force of gravity. There are papers that have been concerned with improving 
the basic Lucas- Washburn model [571 1131 HI]) but here we just use the original model for comparison as a 
reference point. 

As the Lucas- Washburn model predicts that the meniscus is always a spherical cap with, in the case 
of perfect wetting, 9d — 0°, we take the predicted averaged height h, and then fit the spherical cap to it 
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Figure 12: Snapshots of the free-surface shape taken at constant time intervals as predicted by the Lucas-Washburn model 
(left of each figure) compared to our simulation (right of each figure), not to scale. Left: Capillary of radius R, = 0.036cm with 
snapshots at intervals of 100 seconds for the first 2500 seconds. Right: Capillary of radius R = 0.074cm with snapshots at 
intervals of 50 seconds for the first 600 seconds. 

ensuring mass conservation which results in the apex at z a = h— 2R/3 and the contact line at z c — h + R/3. 
This means that at the instant when the simulation starts, the interface jumps from its initially flat shape to 
a spherical cap which extends both above and below the initial interface, i.e. partly back into the reservoir. 
This is a known unphysical feature of the model [58], and forces it to predict infinite speeds as t —> 0. 

In Figure [l2j the predictions of the two models are compared for each capillary. The Lucas- Washburn 
model always predicts a much faster uptake of the liquid, whilst in our simulation the meniscus is slowed as 
it takes time for the dynamic contact angle to decrease, and hence the capillary pressure to increase, from an 
initially flat interface, so that no unphysical jumps in variables are observed in the initial stages. Notably, 
at the end of the time periods considered, the Lucas- Washburn result predicts that the menisci are close to 
their Jurin height, hoo — 1.23 cm and hoo — 0.6 cm, respectively, whilst the results from our simulation are 
such a distance away from that a comparison of the results with experimental ones should be able to 
easily distinguish between the two models' predictions. 

The results from our simulations are compared to both the experimental results of |57j as well as to the 
Lucas- Washburn curve in Figure [13J where two sets of data have been shown in the second graph to illustrate 
the repeatability of the experimental results. Considering that no parameters have been fitted, our compu- 
tations are seen to approximate the data exceptionally well and are vastly superior to the Lucas- Washburn 
curve which hugely overpredicts the speed at which the capillary uptakes the fluid. The inadequacy of the 
Lucas- Washburn equation has been recognized before and attributed to the need to account for a variable 
dynamic contact angle 57J. In a forthcoming article, our code will be used to assess how all the models pro- 
posed for imbibition into a capillary perform, establish the individual effects contributing to the meniscus' 
behaviour and, since none of the models account for the effect of the capillary size on the contact angle, 
ascertain bounds of applicability of the models previously proposed. Here, we arc satisfied with the current 
result which demonstrates the code's ability to describe experimental data with ease. 

In Figure |14[ the relationship between capillary number based on the contact line speed and contact 
angle is presented. Notably, as the two curves are for the same fluid, it is only the variation of the contact- 
line speed that causes the variation of the capillary number, so that the latter can be viewed just as the 



25 



500 



1000 



1500 



2000 



2500 



100 



200 



300 



400 



500 



600 



Figure 13: Apex height of the meniscus (in cm) as a function of time (in seconds). Experimental results |57l are circles and 
crosses (which show a repeated experiment to demonstrate the reproducibility of the data) , our simulations are the solid lines 
whilst the dashed line is the Lucas- Washburn result. 



dimensionless contact-line speed. Surprisingly, the curves for the two different capillary sizes are graphically 
almost indistinguishable. One may naively think that this is an argument in favour of using a speed-angle 
formula; however, as one can see from Figure [M) the velocity dependence of the contact angle is multivalued 
so that it is not a function. For example, at Ca c = 0.03 the contact angle 9d = 45.7° or 84.2°. Furthermore, 
the upper branch of the curve corresponds to a decrease in the contact angle as the contact-line speed 
increases. The upper branch of the curve is far away from the asymptotic line (dashed line) computed from 
(22 1 whilst the lower one is very close to it (the asymptotic value at Ca c = 0.03 is Bd = 46.4°). 

The reason for this double-valuedness of the dynamic contact angle versus the contact-line speed plot 
can be explained very simply. The liquid started moving from rest and at t = the contact angle had 
the value of 8d = 90°. Therefore, if the model does not produce unphysical singularities, with dynamic 
characteristics experiencing instant jumps, then both the contact-line speed and the contact angle have to 
evolve continuously from their initial values, in our case, to and 90°, respectively. Hence there appears 
a branch in the speed-angle plot stemming from the initial state and leading towards the quasi-steady 
state. This manifestation in the speed-angle relationship of the initial conditions rules out in principle the 
possibility of treating the dynamic contact angle as a prescribed input into the model. It is important to 
emphasize that this is a general argument based on the necessity not to have unphysical singularities in 
the mathematical model of a physical phenomenon, and the plot produced by using the interface formation 
model is merely an illustration that a model regarding the dynamic contact angle as part of the solution 
satisfies the above requirement in a natural way. 

Following the initial period, where the dynamic contact angle evolves from its initial value and the 
contact-line speed increases, the flow 'settles' into the 'regular' quasi-steady pattern. There the velocity- 
dependence of the contact angle falls onto the one given by the asymptotics outlined in |2.2| It is important 
to point out here that, although the radii of the capillaries considered in this section are large enough for 
the influence of R on the flow to be negligible, the unsteadiness of the process makes it necessary to use the 
interface formation model in its entirely and not just the asymptotic result for the speed-angle relationship 
that follows from the model in the limiting case. The latter becomes applicable at the later stages of the 



imbibition, whereas at the onset both (22 1 and any of the velocity-dependencies of the contact angle used in 



the conventional models will make the contact-line speed jump from zero straight up to a large finite value. 



7. Concluding remarks 

The necessity to develop a robust scheme of incorporating the model of dynamic wetting as an interface 
formation process into a numerical code is dictated by the need to describe the experimentally observed 
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Figure 14: The relationship between the capillary number based on the contact-line speed (Ca c ) and dynamic contact angle (6^) 
for the two simulations considered, whose curves are graphically indistinguishable, compared to the asymptotic line (dashed) 

effect of the dynamic contact angle being dependent on the entire flow field/geometry. This effect becomes 
more pronounced as the size of the system goes down and has far-reaching implications for a host of emerging 
technologies. 

In this paper, such a scheme is described and the resulting numerical code is thoroughly validated. 
The exposition is purposefully laid out in a 'digestible' manner, with specific details given in the Appendix, 
allowing an interested reader to reproduce the code and compare the results with the benchmark calculations 
provided and, together with our previous publication [57], giving the first guide for simulating fluid flows 
with forming interfaces. A particular emphasis is put on the estimates for the required mesh resolution and 
the physics associated with the smallest length scales for different regimes. These estimates as well as the 
benchmark calculations remain valid independently of the particular numerical scheme used. 

The considered test problem of a steady imbibition into a capillary highlights a unique feature of the in- 
terface formation model, namely that the dynamic contact angle is an output, part of the solution, as opposed 
to being a prescribed function of the contact-line speed and other parameters. For this test problem, our 
computations show that, as the radius R of the capillary is reduced, the velocity-dependence of the dynamic 
contact angle becomes influenced by R, with the dynamic contact angle for the same contact-line speed being 
higher for smaller radii. This is an essentially new physical effect predicted using the interface formation 
model. Given that, when the influence of the flow field/geometry was first discovered experimentally |15j . 
the term 'hydrodynamic assist of dynamic wetting' was coined to emphasize the possibility of lowering the 
dynamic contact angle by manipulating the flow field, here we may term the newly discovered phenomenon 
as 'hydrodynamic resist to dynamic wetting'. The focus of this paper is on numerical implementation of 
the interface formation model and so we postpone the full exploration of the effect of 'hydrodynamic resist' 
to a future publication in which our framework will be compared to experiments, e.g. [17] . where resist-like 
behaviour has been observed but is yet to be theoretically described. 

The problem of unsteady imbibition into a capillary demonstrates the above unique feature of the inter- 
face formation model from a different perspective. This problem shows that the velocity-dependence of the 
dynamic contact angle bears the influence of the initial conditions even in the situation where otherwise the 
effect of the flow field/geometry would be negligible. In particular, the dynamic contact angle can decrease 
with the increasing contact-line speed. Qualitatively, this possibility follows from a general argument about 
the continuity of variation of the flow characteristics necessary to avoid unphysicality, and the interface 
formation model provides a quantitative illustration of how this unphysicality is removed. 

The validation of our numerical code against experimental data demonstrating excellent agreement with- 
out the need to fit any parameters opens up two lines of inquiry. First, it would be interesting to consider 
how scaling down of the system towards submicron and nanoscales affects the flow and, by comparing the 
results with experimental data, identify the scales where specific 'nanoeffects' become important. Secondly, 
the code can be used to validate different simplified models for the capillary imbibition by identifying the 
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regions of parameters where different effects play the dominant role. 

What has been shown is only a small part of the capabilities of the interface formation model fully 
incorporated into a numerical framework. In [29], the dynamics of freely oscillating liquid drops and of 
drops impacting and spreading on chemically patterned surface is considered, demonstrating the framework's 
capability to simulate flows in which changes in free-surface shape are extreme. The ability of the interface 
formation model to incorporate additional physical effects is shown in [371 160) , where a viscous flow over a 
solid surface whose wettability varies has been described, with the variations in wettability affecting the flow 
even in the absence of any meniscus. Further influences on the wetting dynamics, such as those associated 
with the solid substrate's properties, e.g. its porosity [ST], or thermal and electromagnetic effects, can be 
added as required, safe in the knowledge that the underlying framework's accuracy has been completely 
confirmed. 



8. Appendix 

In |27j , a user- friendly step-by-step guide to the finite element simulation of dynamic wetting flows using 
the conventional models is given. Many of the aspects dealt with in that paper remain unchanged when we 
introduce the interface formation model into this framework, and these are not repeated here. In particular, 
we do not repeat: the re- writing of tensorial expressions in terms of a particular coordinate system, the mesh 
design which utilizes the bipolar coordinate system, the mapping of integrals onto a master element, the 
calculation of the corresponding integrals using Gaussian quadrature or the solution procedure based on the 
Newton-Raphson method. Here, the main change is that we are now solving many more equations along the 
boundaries of the domain, and, consequently, the element-level residuals change from those in [27] and are 
thus listed in full. Therefore, in parallel with [37], this guide enables one to easily reproduce all the results 
presented in the paper. For a more detailed exposition about using the FEM to solve fluid flows the reader 
is referred to 62, 63 , whilst specific information on free surface flows, including alternative mesh designs, 
can be found in [64, 65, 66, 67, 50]. We have provided information of our own specific implementation of the 
framework described in the main body of the paper. Many alternatives exist (e.g. different element types) 
which may produce equally accurate results but here we present what we have used and tested, which for a 
'practitioner' gives a ready-to-use algorithm. 

Two-dimensional and three-dimensional axisymmctric flows are considered simultaneously by using a 
variable n which takes the values n = in the former case and n — I in the latter. Two-dimensional flow 
occurs in a Cartesian coordinate system (r, z) whilst axisymmetric flow is in a cylindrical polar coordinate 
system [r, z, 1?), where r is now the radial coordinate and $ is the azimuthal angular coordinate around 
the z-axis about which symmetry is assumed (e.g. Figure [2]). For both coordinate systems, the governing 
equations are solved in the (r, z)-plane. In this coordinate system, the surface velocity is v s = vft + v^n, 
and vfj = uf t as, in the n = 1 case, vfj ■ eg = 0. Then, notably, the divergence of the tangential component 
of a surface vector is given by 

yjs » . g (aft) , n < da} na s r 

V ■a,i=t-— 1 = — H (53) 

11 os r os r 

dt 

where we have used that t • — — = and at — & s ■ e r . 

os 

8.1. Finite element procedure 

The finite element method presented here is not restricted to the capillary geometry and our only 
assumption about the free surface is that in the (r, z)-plane, i.e. the computational domain, it is a line 
parameterized by arclength s with position defined by a function of one variable h = h(s) which depends 
on the particular mesh design. In other words, to determine the free-surface shape we must determine the 
value of hi at every node on the free surface i = 1, Nx . 

The domain is tesselated into e = 1, N e non-overlapping subdomains of area E e , the elements, which 
in our case will be two-dimensional curved-sided triangles whose positions are defined by the spines of the 



28 



mesh and hence by the values of free surface unknowns hi. The boundary of the computational domain is 
composed of one-dimensional elements of arclength s e formed from the sides of the bulk elements which are 
adjacent to the boundary. We refer to element-level quantities as local and those defined over the entire 
domain, which were used in Sj3j as global. Each element contains a set number of local nodes, and it is the 
set of all local nodes which form the global nodes referred to in f|3] We have used Roman letters (i,j) for 
global quantities and italicized letters (i,j,k,l) for local ones. In the FEM, one must store a function I 
which relates each local node i in each element e to its global node number i so that i = J(e, i). 

The global functions and residuals are constructed in a piecewise manner from local quantities, by 
ensuring that, as required in the FEM, the interpolating function associated with a given node is constructed 
to be zero outside the elements to which that node belongs. This allows us to derive expressions for the 
residuals in an arbitrary element which, after summing up contributions from every element in the domain, 
will form a set of algebraic equations whose coefficients are integrals to be calculated by Gaussian integration 
as shown in [27] , 

To avoid having to construct each local interpolating function in an element-specific manner, they are 
constructed on a master element with coordinates (£, 77), so that <pi — (f>i(£,T]), ipi — V , i(Ci 7 ?)- Mixed inter- 
polation is achieved, to ensure the Ladyzhenskaya-Babuska-Brezzi [45 condition is satisfied, by using the 
V6P3 Taylor-Hood triangular element which approximates velocity by means of bi-quadratic local interpo- 
lating functions 4>i(£,,r]) (i — 1 , • ■ ■ , 6) and pressure using bi-linear ones ipi(^,r/) (i = 1,...,3), with explicit 
expressions given in |27j . Over an element e the following functional forms are used 

6 3 6 

(u,w) =^T(u j ,w j )(t> j (£,Ti), P = ^Pj*l>j(€,V)> ( r > z ) = ^2( r 3> Z j) ( i ) j(€>V)- 
i=i i=i j=i 

By ensuring that elements whose sides form the free surface always have local nodes 2, 6, 3 associated with 
them and that those on the solid surface are always j = 1, 5, 2, we can define surface interpolating functions 
(f'l.jiv) = 'Pji^V = — 1) an d 4>2,j(j() = 4>j(£ = — 1 , ??) , respectively. Then, a surface variable a s and the 
free-surface shape (ri(h), Zi(h)), which is determined from the free surface unknowns h, are approximated 
on the surfaces quadratically as 

J=2,6,3 i=2,6,3 J=l,5,2 

The global residuals (the Ri's) presented in ^3j which involve integrals over the entire domain, are formed 
by summing up local residuals (the R e /s), denoted with a subscript e to indicate which element the local 
residual is calculated in, obtained by integrating over each of the N e elements in the global domain. Let 
Ni be the set of numbers of those elements that form part of the free surface and N2 be the numbers of 
elements forming a part of the liquid-solid interface. Then, the continuity of mass, momentum equation^ 
kinematic equation on the free surface and impermeability equation at the liquid-solid interface are: 

N G 3 N e 6 W e N e 

e=l i=l e=l i=l e=l i=2,6,3 e=l i=l,5,2 

i=/ ( e ^) i=I(e,i) i=I(e,i), eSJVi i=I(e,i), e£W 2 

where the constraint under the summation symbols ensures that, through i = I(e,i), the local residuals are 
assigned to the correct global ones and, via e € N, that the surface equations are only evaluated when an 
element is adjacent to the boundary of the domain. In addition, we now have residuals from the interface 



7 Note the typo in |27| where the limit for i in Rp and R^' a should be 3 and 6, respectively 
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formation equations given by 



R 



R 



E E 

e=l i=2,6,3 
i=7(e,i), c£Wi 



E E 

e=l i=l,5,2 
i=J(e,i), eG7V 2 



JV, 



EE *£f. *T 

e=l i=2,6,3 
i=I(e,i), eeNi 

N e 

E E *Sr. *f' 

e=l i=l,5,2 
i=J(e,i), eeS 2 



E E R r 

e=l i=2,6,3 
i=J(e,i), eeJVi 



E E 

e=l i=2,6,3 
i=I(e,i), eeJVi 



P 2 



7V e 

E E 

e=l i=l,5,2 
i=J(e,i), e£N 2 



R P2 - R a 



(55) 

E E *3> 



A', 



e=l i=l,5,2 
i=J(e,i), eGJV 2 



(56) 



and the residual from Young's equation R Y which is only applied at the contact line node. 

Taking the equations of ^j3]over an arbitrary element e, with area E e and boundary s e , we now derive 
the local residuals required in (|54| and ([55l. 



8.2. Element-level residuals 

For an arbitrary element e, with local nodes i = 1, 
to the momentum residuals are 



.,6, contributions from the bulk of the clement E e 



R™f = Rc 



R™f = Rz 



3 dt 

3 dt 



■A, 



A,, 



dt 



dt 



Kjf Wj 



Kffa + K%w 3 + Cl Pk + St d, (57) 
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1. 



,6 k = l, 



where summation over repeated indices is henceforth assumed. The M terms are the mass matrices, the A 
terms are from the nonlinear convective terms, K terms are associated with viscous forces, the C terms are 
with pressure forces and G terms correspond to the body force acting on the liquid due to gravity. Temporal 



derivatives will be considered in £8.3 



For i = 1, ...,3 we have the incompressibility residuals 

Re.i — CjiUj + CjiWj 



3 = 1,. ..,6. 



If an element e forms a part of the free surface, e € N\, then for i = 2, 6, 3, i.e. the free surface nodes, 
there are additional contributions to the momentum equations from capillary stress terms, so that 



M,r 



A/.; 



M,r 
e,i 
M,z 
e,z 



h F ij a ^ ~ In \jPs,o 3 = 2 > 6 > 3 . 



Ca y 

m F ij a ^3 ~ In ijPgJ 



3 =2,6,3, 



(58) 



where p g j is the nodal value of the gas pressure, which in this work has been assumed to remain constant. 
If an element e forms the part of the free surface, e G JVi, adjacent to the contact line, so that local node 
i = 2 is the contact-line node, then there is an additional contribution 



H e,2 ~ K e,2 



R 



M,r 
e,2 

M,z 
e,2 



--R. 



M,z 
e,2 



1 T 1 ^ 

C^ 1 °V=2: 

1 ^2 



(59) 



which is where the contact angle is imposed into the weak formulation. The contact angle is determined 
from the Young equation, now applied only at the contact-line node, whose residual is given by 



i? 1 



ai,i=2 cos6> d + <7 2 ,i=2- 
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Additionally, for i = 2, 6, 3 one has the kinematic equation 



R. 



K 

e.i 



In 



dt 



dt 



J =2,6,3. 



If an element e forms a part of the liquid-solid interface, e € A^2, then for £ = 1,5,2 there are additional 
contributions to the momentum equations of 

Rff =R J ff + JL ( N £( Uj - U : ) + Ngfa - Wj)) + In^Aj - ^Sfaj j = 1, 5, 2, 



R M f =R 



e,i ' Ca "11 ' "]') 1 2Ca V 

M . z i / i\r21/„, rr ^ i (\r22/„,, tj/ ^ i r„2 a 1 o2 



+ A - ^) + ^ K - ^)) + J 4 A J - J = 1. 5- 2, (60) 



where the N terms come from the slip terms in the generalized Navier condition ( 33 1 , whilst the 5 terms are 



associated with gradients in surface tension in ( 33 ) and the terms containing A are from the normal stress 
on the interface. The substrate speecj^] has been decomposed into scalar components as U = Ue r + We z . 
Additionally, for i = 1,5,2 one has the impermeability equation (46) in the form 



R\, = kivlnj ~ {HjUj + InlWj) j = 1, 5, 2. 
Moving on to the equations of the interface formation model, which do not appear in [27] . the Darcy-type 



equations ( 36 1 and ( 37 1 give 



r: 



(lt\ jUj + It%w 3 ) - > *, J = 2, 6, 3, 



R V e,i = - 3 [ R U U J + U j) + n U w 3 + W i)] ~ *Dii°2J hi = 1» 5 > 2- 

The surface velocity normal to each interface 7 = 1, 2 is determined from 

Kj = ( In ij u i + In % w j) ~ tijV^nj - Qlij (p 7 j ~ P s 7e ,j) 7=1: i,j = 2, 6, 3, 7 = 2: i,j =1,5,2 
where we have allowed for the extension to problems in which varies spatially. 



The terms appearing in the surface mass balance equation ( 41 1 are more complicated due to the presence 



in this equation of temporal derivatives and the divergence of surface vectors. First, recall that on each 
surface 7 = 1,2 the finite element surface mass balance equation is 



R 



7 |M ' VV., 



m 7 dC 7 



-gf + PjV- C 7 || + P7K ' n 7) V ' n 7 ) + ^74 {(£, ~ Pye) 

(i = l,...,A 7 ). (61) 



Before deriving the individual terms, it proves convenient to introduce the component of the surface's mesh 
velocity tangential to that surface as a new variable which we label c s t . This is discretized in the usual way 
on each surface 7 = 1,2 



7: 



This new variable is in the (r, z)-plane, so that cf 

dz 



R 



7* 

e.i 



T r s 

y 7*J 



1 dr,- 
13 dt 



dt 



1: j = 2,6,3, 7 = 2: j = 1,5,2. 

= cfi • t and is thus defined by the finite element equation 
7 = 1: i,j =2,6,3, 7 = 2: i,j = 1,5,2. 



Which should also have been present in 1271 . 
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The first of the three nonlinear terms in (61) is 
ds 



"'V" ( V 'jt ~ C 7t) rn ds ie = P^j (« 7 t,fc ~ C*t tk ) 



d<p. 



PijkP^j ( v yt,k c 7*,fe) — Pij ( v yt,k> c jt,k) P7 



(62) 



The second one is 



7 ' fe r"ds 7e + np. s dri ' k 



W 7J ds 

w t3k pt ( jc: (tik + nw- jk p^ tj 



™ di 
dn. 



dt 



^7, i^7,i ( ^ > 7, A; ds~^ e 

WiiC^fc,^,*)^. (63) 



The third term gives 

/ ^^^^V 8 • n 7 r" ds 7e = p 7J w 7 „, fe / <t> 1 ,i4> 1 ,A 1 .kKr n ds -ye = VijkPyjV^ = Y^iv^^p^^ 



dn ir 

where k — t 



■U 



di 



(*S). 



1 n(n ir /r) is the curvature, 
lass 



7I " 9s 9s "V'fcyr 

Finally, this gives the surface mass balance equation as 



«7 



dt 



I'uP-.j + W ijPj,J + 



Whenever the interface reaches the contact line, where now mi = ti and r = r c , we have an additional 
contribution on both the free surface and solid surface sides as: 



(< : ) C „ = -W 



V lt,i=2 C lt,i=2) r c J c [ 



Lastly, the equations of state ( 46 ) are given by 



where h = Y^j ^ 



8.2.1. Element-level matrices 

It should be pointed out that, in contrast to [27], in this paper non-dimensional parameters are not 
contained in the element-level matrices below as they have already been factored out. 

The mass matrix is 



M„ 



and temporal derivatives are handled in £8.3 
The viscous terms are 



K\] = 2fc.„n + fc i?2 2 + 2n fcj*-, Kif = %2i, Kfj = k^u, Kff = L- u + 2k 



-12 



r n dE„. 



21 



^22 



where 



^ijkl 



%"22 7 



r" dE e 



_ 9t/ fc % 

with j/i = r, j/2 = z. Nonlinear convective terms are 

AjKH = a vj kiu k + a ljk2 w k k = 1,...,6, 
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a i3 ki= I fafik^r 1 - r n dE K . 



where 

n .* „■ j. i — / rh:fh>. 

dyi 

The C matrices represent pressure terms whilst their transpose is required in the continuity of mass equation: 

The gravitational force contributes through the G term and is given by 

Gi= [ (f> t r n dE e . 

On a free surface (one-dimensional) element s\ e , the capillary stress terms are 

rpl [ . (d<t>i,i . nfaA 2 f d ^,i, n, 

F ij = J *lr y-^T + j <Kj r ds le, Fy = J tu-j^-faj r ds U . 

The components of the vectors tangential t = (t^ r ,t 7Z ) and normal n = (n ir ,n lz ) to surface 7 = 1,2 are 
found in the usual way, see [27] . 

At the contact line (r c , z c ) contributions will occur at local node i — 2 of the first free surface element 
so that 

T 1 = [(t 2r cos 9 d + n 2r sin d ) r n c ] , T 2 = [(t 2z cos 9 d + n 2z sin d ) r™} . 

On a (one-dimensional) clement s 2e on the liquid-solid interface, the tangential stress terms from the 
Navier-slip boundary condition for k,l — 1,2 are 

N ij = h,ih,jhkt 2 ir n ds 2e , 

with t 2 i = t 2 r and t 22 — t 2z . The terms associated with gradients in surface tension are given by 

J <P 2 /-^t 2r r n ds 2e , &,,^t2,r»d* 3ej 

On a surface 7 = 1,2 the following other expressions have been used 



/ 7 ,i^ 7 ,j^ 7 r ^* ds*y e7 ^^ij / ^ 7l 'i0 7 ,j^ 7 2 ^* ds^ e . 



where 



Dij — j <j> lt i r ds ie , Pij — Pijk ( v ^t,k c 7 t,fc) 1 

Pijk — I g s 7 ,j<p 7 ,fc r ds Je . 

nl _ I U V-<A 1 1 . n j q 2 _ [ @&r>i , , . n j 
^ij'A: — y v ? 7 ,j ( ?7,fcf' 7 r 71 "S 7e , — y ^7 ,j '^P^.k^fz T dS^ e 

tut s r dr 7 £■ 

Wij = w l]k c ltk + nw ljk -^-, 
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where 
where 

Uijk = / <l>-(,i4>-y,j4>-y,kKr n ds 7e . 



8.3. Temporal discretization 

The result of our spatial discretization is a system of non-linear Differential Algebraic Equations (DAEs) 
of index 2 [63]. However, it is well known that the same methods that apply to ODEs can be used for DAEs 
[68] . and for a review of the methods available we refer the reader to [53]. 

The second-order Backward Differentiation Formula (BDF2), which has been applied successfully to 
similar problems [41] . is implemented into our scheme. Below, the method is applied to a scalar equation 
y = f(y,t), with the extension to derivatives in the Navier-Stokes and interface formation equations being 
a straightforward task. For a time (n + 1) with step At, the method applied to the scalar equation gives 

Vn+l - Vn _ \ Vn~ Vn-1 2 . 

At ~ 3 At + 3 y " +1 ' 

where the subscript indicates the time step at which a variable is evaluated and dy/dt = y. Alternatively, it 
may be written as 

3y n +i - + yn-i _ . 

2At -Vn+U 

which is the second-order accurate one-sided Taylor series expansion of y at t n +i- 

Often, during a physical process there will be different stages that are characterized by different time 
scales. It is important that our temporal discretization takes these different scales into account so that the 
largest possible time step, that enforces a certain accuracy, is chosen automatically. This is achieved by 
choosing a 'step' so that the local truncation error d n = y n +i — y(t n +i), where y(t n+ i) is the exact solution, 
is maintained below a certain tolerance. By using an explicit second-order Adams-Bashforth method (AB2) 
[62] to predict the solution y n+11 and then comparing the difference between the actual solution and the 
predicted one, after solving the non-linear equations we are able to deduce this error. For the equation, 
y = f(y,t), having obtained the solution y n +i we follow the analysis of |62j to obtain a new time step. The 
BDF2 for variable step size gives 

Vn+l - Vn _ At n y n - y n -i At n + Atn-l . . . 

At n ~2Ai„ + At„_! At„_ x 2At n + At n - 1 yn+1 ' [ ' 

while the predictor gave 

At, n \ . , > ( At 



y'n+i = Vn + ( 1+ At** 1 1 ) ~ ^™ ~ Vn ~^ 



At„. 



The local truncation error of ( 64 ) is given by 

(At n 

A*„(2Ai n + At n -i) 6 



(At n + At n ^) 2 Atjy n ^ fA ^ 

d n = a 4- /o A —7T+ \ a + °( A ^) (65) 



whilst the predictor's error is 

„ p . . - »/f*_ ,,1=-IH 

At n -! J 6 



y p n+1 - y(t n+ i) = ( i + ~) ^ + o(A<). (66) 
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The exact solution's contribution y(t n +i), which, of course, in general will not be known, is eliminated from 
( 65 1 and ( 66 ) to give 



, = (1 + At w _i/At n ) 2 _ p 

" 1 + 3(At n _i/At n ) + 4(At„_i/At n )2 + 2(Ai n _!/Ai n )3 {Vn+1 Vn+lh 

so that the error is linearly proportional to the difference between the predicted and actual solution. This 
error estimate is then used to compute the next step size 

At n+1 = At n (e/\\d n \\) 1 / 3 , K|| 2 = dld n /{Ny 2 max ). 

Here, N is the total number of nodes, y m ax is an estimate of the maximum value of y in the domain 
and e is the relative error tolerance parameter. Then the error of the approximate solution is bounded by 
ll^n+ill < tymax- This allows us to choose the largest possible time step whilst ensuring that the error of 
the temporal integration remains below a chosen tolerance. 

The variable step method outlined above can be extended for use with the Navier-Stokes and interface 
formation equations and the reader is referred to [631 p. 797], for details. In this book a number of 'rules 
of thumb' are suggested based on the ratio of the new and previous step sizes, the so-called Delta T Scale 



Factor, DTSF = 



and these can be adapted to the problem of interest as required. 



8.4- Far-field velocity profile for the steady propagation of a meniscus through a capillary 

To ensure mass is conserved in our domain for this steady problem in which an adsorption-desorption 
process will occur along the interfaces, we must derive appropriate conditions at the base of the capillary, 

dp s 

i.e. in the truncated far field. Assuming the problem is steady, so that Cy = = 0, noting that there is 

no surface mass flux through the axis of symmetry and that there is no sink or source of mass at the contact 
line, we have 



= 



/ V • udV = / u- ndA = / u ndA ff + Q {p\ - p{ e ) dA x + Q (p| - p s 2( 

JV J A JAf f J Ai JA 2 



) dA 2 



wdA 



(27T)" 



wr n dr — eQ 



s s 

Pl v ll 



m.\r 



apex 



ff 



Pi v i|| • mir 




with subscripts // referring to far field variables. 

Given that w = w(r) in the far-field, from the equations of motion we see that 



dp 
d^ 



1 d 

r n dr 



, dw 



and, after integrating once, with 



dp 

dz 



G, this gives 
Gr 



Ar- n = 



Noting that 



dr 



1 ' " dr' 
at r = 0, so that A = 0, we integrate again to find 

&3 + B. 



2(n+l) 
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Applying the Navier condition — — = (f3/Ca) (w — W) at r = 1, assuming that in the far field there are 

dr 

no gradients in surface tension, we have that 

G E 



n + 1 Ca \2(n+l) 



G 



so that 



and, finally 



B = W- 



w = W + 



G 

2(n + l) 
G 



2(n + l) 



r 2 -l 



/3/Ca 
2 



W/Ca) 



Now, we calculate the pressure gradient, G, required to maintain a steady flux in and out of the domain 
for a given W by calculating 



wr n dr 



„n+l 



W- 



G 



2r 



n+1 



n + 1 2(n + l)\n + 3 n + 1 (f3/Ca){n + l) 



-i l 



r=0 



w 



1 



G 



n + 1 2(n + 3)(n+l) 5 



n + 1 — n — 3 — 



2(n + 3) 
ifi/Ca) 



(67) 



Noting that the surface velocity in the far field is v^t = (1/2) (w + W), where q — £Qp?, e , we have that 



= W 



1 



G 



n+1 (n + 3)(n + l) s 



n + 3 



ifi/Ca) 



qW 



Gq 



2{n + l)(l3/Ca) 



so that 



G 



2W{n + l){l + q) 
2/(n + 3) + {2 + q)/(/3/Ca) 



with q = q(n + 1) and therefore finally 



w = W + 



W{1 + q) 



2/(n + 3) + (2 + q)/{0/Ca) 
where, in f|4j we had n = 1 and W = — 1. 



(r 2 -l-2/(^/Ca)) 
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